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Abstract 

The  overall  objective  of  this  research  effort  was  to  formulate  a  preventive 
maintenance  strategy  for  AMRAAM  missiles  subject  to  extended  captive  carry  flight 
time.  A  preventive  maintenance  policy  is  only  applicable  if  the  item  in  question  is  aging, 
or  deteriorating  with  time.  Therefore,  a  supporting  objective  of  this  research  is  to 
characterize  the  aging  process  of  the  missile  system  through  a  non-parametric  analysis  of 
its  Mean  Residual  Life  (MRL)  function.  Three  non-parametric,  censored-data  MRL 
function  estimation  techniques  discussed  in  the  literature  are  examined  via  a  numerical 
example.  All  three  estimation  techniques  provide  MRL  functions  that  exhibit  greatly 
exaggerated  decreasing  trends  compared  to  the  MRL  function  of  the  underlying 
distribution  in  the  example.  A  semi-parametric  technique  for  estimating  the  MRL 
function  is  developed  that  shows  dramatic  improvement  over  the  non-parametric  results. 
Although  the  MRL  analysis  of  the  current  AMRAAM  failure  data  failed  to  provide 
evidence  that  the  missile  system  is  aging,  three  preventive  maintenance  policies  discussed 
in  the  literature  are  investigated.  The  traditional  approach  of  preventive  maintenance 
policy  optimization  via  cost  function  minimization  requires  the  cost  of  a  system  failure  be 
explicitly  known.  However,  the  penalty  for  a  system  failure  is  often  subjective  and 
difficult  to  express  in  monetary  terms.  A  “reliability  cost  model”  is  developed  whereby 
system  reliability  for  each  policy  is  expressed  as  a  function  of  cost.  This  technique 
allows  a  direct  assessment  of  the  trade-off  between  cost  and  projected  system  reliability. 
Theoretical  results  are  presented  and  the  performance  of  the  model  applied  to  empirical 
data  is  assessed  via  a  Monte  Carlo  simulation. 
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OPTIMUM  PREVENTIVE  MAINTENANCE  POLICIES 
FOR  THE  AMRAAM  MISSILE 


1.  Introduction 

Overview 

The  Advanced  Medium  Range  Air-to-Air  Missile  (AMRAAM)  is  a  highly  capable 
active  radar  guided  missile  carried  by  USAF  F-15  and  F-16  fighters.  The  AMRAAM 
missile  system  is  composed  of  four  major  components:  1)  the  seeker  section,  2)  the 
warhead  /  fuse  section,  3)  the  flight  control  section,  and  4)  the  rocket  motor  section.  The 
missile  has  a  built-in-test  (BIT)  capability  whereby  the  major  missile  components  are 
electronically  monitored  and  failures  are  reported  to  the  pilot.  When  a  failure  is  indicated 
the  missile  is  no  longer  launch  capable. 

Operational  requirements  over  the  last  several  years  have  caused  the  missiles  to 
accumulate  extensive  “captive  carry”  hours  whereby  the  missiles  are  carried  on  the 
fighters  without  being  launched.  The  AMRAAM  was  initially  designed  for  a  captive 
carry  mean  time  between  failure  (MTBF)  of  1000  hours.  However,  a  low  frequency 
vibration  problem  discovered  when  the  AMRAAM  is  captive  carried  on  fuselage  stations 
of  the  F-15  prompted  the  government  to  reduce  the  MTBF  acceptance  requirement  to  450 
hours.  Due  to  current  mission  requirements  most  AMRAAMs  being  flown  overseas  have 
in  excess  of  450  hours  and  many  have  in  excess  of  1000  hours  captive  carry  time  and  a 
large  number  of  BIT  indicated  failures  during  captive  carry  operations  have  been 
recorded. 


1 


Williams  and  Pohl  (1997)  analyzed  AMRAAM  failure  data  and  determined  the  life 
time  distribution  of  the  missile  does  not  always  exhibit  a  constant  failure  rate.  In  fact, 
several  periods  of  increasing  failure  rate  were  found  indicating  the  missiles  age,  or  wear, 
with  respect  to  cumulative  captive  carry  time.  However,  the  exact  number  of  captive 
carry  hours  at  which  the  missile’s  performance  is  degraded  remains  an  open  issue  and  no 
preventive  maintenance  policy  is  currently  being  used  to  improve  performance. 

Background 

The  objective  of  a  preventive  maintenance  policy  is  to  minimize  system  failures  via 
replacement  of  components  that  reach  a  life  at  which  the  system’s  reliability  gets  to  be 
lower  than  the  reliability  goal  set  for  the  next  mission.  (Kececioglu,  1995:  xxviii) 
Reliability  is  defined  as  “the  probability  an  item  will  adequately  perform  its  specified 
purpose  for  a  specified  period  of  time  under  specified  environmental  conditions.” 
(Leemis,  1995:  2)  Two  criteria  must  be  met  before  a  preventive  maintenance  policy  is 
considered.  First,  the  item  in  consideration  must  be  aging,  or  deteriorating  with  time.  If 
the  item  remains  “as  good  as  new”  regardless  of  age  (exponential  lifetime),  or  is  “used 
better  than  new”  (improving  with  time),  then  preventive  maintenance  will  not  improve 
reliability.  Second  the  cost  of  a  system  failure  must  be  greater  than  the  cost  of 
performing  preventive  maintenance.  If  the  penalty  of  an  in  operation  failure  is  no  greater 
than  the  penalty  for  performing  scheduled  maintenance  then  no  cost  benefit  is  realized 
even  if  reliability  is  improved.  Given  these  criteria  are  met,  the  best  preventive 
maintenance  policy  is  one  that  optimizes  some  objective  such  as  minimizing  cost, 
maximizing  availability,  or  maximizing  reliability. 


2 


The  reliability  function,  S(t)  (sometimes  called  the  “survivor  function”)  of  an  item  is 
given  by  S(t)  =  1  -  F(t)  where  F(t)  is  the  lifetime  cumulative  distribution  function.  When 
the  item  of  interest  is  a  system  composed  of  several  independent  components  such  that 
the  failure  of  any  one  component  causes  a  system  failure  (a  series  system),  the  system 
reliability,  So{t)  can  be  expressed  in  terms  of  the  component  reliabilities: 

nc 

■So(')=n*.«  (1) 

1=1 

where  Si(t)  =  reliability  of  component  f 

nc  =  number  of  components. 

Typical  maintenance  models  require  specific  knowledge  of  the  survivor  functions  at  both 
the  system  and  component  level.  Estimates  of  the  survivor  functions  must  be  made  from 
available  failure  data  that  is  often  censored.  Data  censoring  occurs  when  specific  failure 
times  are  not  known  for  all  items  on  test.  Right  censoring  occurs  when  an  item  is 
removed  from  test  before  failure  for  some  reason.  In  this  case  only  a  lower  bound  is 
known  for  the  item’s  lifetime.  There  are  three  types  of  right  censoring.  Type  I  censoring 
occurs  when  the  test  is  terminated  at  some  specified  time,  to-  The  failure  times  for  all 
items  that  have  not  failed  by  this  time  are  censored  with  a  lower  bound  of  to-  Type  II 
censoring  occurs  when  the  test  is  terminated  at  a  specified  number  of  failures.  Type  III, 
or  random  censoring,  occurs  when  items  are  removed  from  test  at  any  time.  The  current 
AMRAAM  failure  data  set  consists  of  a  large  number  of  observations  (n  =  815)  subject  to 
approximately  74%  random  right  censoring.  Future  data  is  projected  to  have  over  2500 
observations  subject  to  approximately  80%  censoring. 
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Research  Objectives 

The  overall  objective  of  this  research  effort  is  to  formulate  a  preventive  maintenance 
strategy  for  using,  retiring,  or  refurbishing  AMRAAM  missiles  subject  to  extended 
captive-carry  flight  time.  This  objective  requires  assessment  of  the  two  aforementioned 
criteria  for  implementing  a  preventive  maintenance  policy.  Namely,  is  the  system  aging 
and,  if  so,  is  the  cost  of  an  in-operation  failure  greater  than  the  cost  of  performing 
preventive  maintenance?  We  assume  the  answer  to  the  second  question  is  yes,  even  if  the 
penalty  for  an  in  operation  failure  is  not  measured  directly  in  dollar  terms.  The  aging 
criterion  is  yet  to  be  resolved.  Williams  and  Pohl  (1997)  found  evidence  of  aging 
through  periods  of  increased  hazard  rate.  However,  the  exponential  distribution,  which 
has  a  constant  failure  rate  and  is  indicative  of  an  item  that  does  not  age,  was  a  good  fit  to 
the  failure  data.  Therefore,  a  supporting  objective  of  this  research  is  to  characterize  the 
aging  process  of  the  missile  system  through  a  non-parametric  analysis  of  its  Mean 
Residual  Life  (MRL)  function. 

Scope 

We  make  a  simplifying  assumption  that  there  are  infinite  spares  available.  In  other 
words,  if  a  missile  is  taken  out  of  service  for  maintenance,  planned  or  otherwise,  a  spare 
missile  is  always  available  to  fill  mission  requirements.  Under  this  assumption 
preventive  maintenance  does  not  affect  missile  availability  and  the  time  required  to 
perform  preventive  maintenance  need  not  be  considered  in  the  cost  models. 

Overview  of  Subsequent  Chapters 

Chapter  2  contains  a  review  of  the  existing  literature  covering  the  topics  pertinent  to 
this  research.  Non-parametric  estimation  of  the  MRL  function  and  tests  for  aging  are 
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explored.  Several  preventive  maintenance  policies  with  their  associated  cost  models  and 
optimization  techniques  are  also  examined.  The  methodology  for  the  final  analysis  is 
developed  in  Chapter  3.  Existing  techniques  for  MRL  analysis  are  validated  and  a 
technique  for  improving  the  non-parametric  estimate  of  the  MRL  function  with  censored 
data  over  those  described  in  the  literature  is  developed.  A  nontraditional  approach  to 
maintenance  policy  optimization  is  developed  whereby  the  objective  is  to  maximize 
system  reliability  subject  to  an  unspecified  cost  constraint.  The  methodologies  developed 
in  Chapter  3  are  applied  to  the  AMRAAM  failure  data  set  and  results  are  reported  in 
Chapter  4.  Chapter  5  contains  a  summary  of  the  thesis  effort,  including  an  overview  and 
discussion  of  the  impact  of  the  results  and  ideas  for  future  research. 
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11.  Literature  Review 


Overview 

This  chapter  provides  an  overview  of  the  current  literature  in  areas  pertaining  to  this 
thesis.  The  chapter  begins  with  a  review  of  the  MRL  function  and  its  applications.  This 
review  includes  aging  analysis,  non-parametric  estimation  of  the  MRL  function,  and 
statistical  tests  for  the  decreasing  mean  residual  life  (DMRL)  and  new-better-than-used  in 
expectation  (NBUE)  life  distribution  classes.  The  next  section  demonstrates  the  various 
MRL  concepts  via  a  numerical  example  that  includes  both  complete  and  censored  data 
cases.  In  the  final  section,  three  common  preventive  maintenance  policies  and  their 
associated  cost  models  are  examined.  This  section  concludes  with  a  discussion  of 
maintenance  policy  optimization. 

Mean  Residual  Life 

The  MRL  Function.  The  mean  residual  life  (MRL)  function,  m{t),  is  defined  as  the 
expected  remaining  lifetime  of  an  item  given  the  item  has  survived  to  time  t.  Let  Tbe  a 
random  variable  denoting  the  lifetime  of  an  item  with  associated  cumulative  distribution 
function  (cdf)  F(t),  survivor  function  S{t)  =  \-F(t),  density  function  fit),  and  hazard 
function  h{t)  -  f(t)fS(t) .  By  definition,  the  MRL  function  can  then  be  expressed  as: 

mit)  =  E[T-t\T>t].  (2) 

Like  F(t),fit),  and  h(t),  the  MRL  function  is  unique  to,  and  completely  determines,  the 
distribution  of  T.  The  MRL  function  can  be  expressed  in  terms  of  the  other  lifetime 
distribution  representations  and  vice  versa.  A  common  representation  of  the  MRL 
function  is  in  terms  of  the  survivor  function  as  shown  by  equation  (3). 
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(3) 


fn(t)  =  -^]siu)du. 


Relationships  between  the  MRL  function  and  other  lifetime  distribution  representations 
are  shown  in  Table  1. 


Table  1.  Lifetime  Distribution  Representations  (Leemis,  1995:55) 
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Note  that  m(0)  =  J  uf(u)du  is  the  expected  value  of  T. 

0 

The  Weibull  distribution  is  commonly  used  in  lifetime  modeling  due  to  the  flexibility 
afforded  by  its  two  parameters.  The  shape  parameter  /3  controls  the  shape  of  the 
distribution  and  allows  for  aging,  memoryless  (exponential),  or  improving  properties. 
The  scale  parameter  a  controls  the  spread  of  the  distribution  across  the  time  axis. 
Weibull  distribution  representations  (Leemis,  1995:88)  are  shown  in  Table  2,  and 
Weibull  survival,  density,  hazard,  and  MRL  functions  with  various  shape  parameters  and 
scale  parameter  a  =  1  are  plotted  in  Figure  1. 
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Table  2.  Weibull  Distribution  Representations 


A  Weibull  distribution  with  shape  parameter  =  1  reduces  to  the  exponential 
distribution.  Note  that  in  this  case  both  the  hazard  function  and  the  MRL  function  are 
constant.  Shape  parameter  jS  >  1  indicates  a  deteriorating  item.  Note  the  hazard  function 
is  increasing  and  the  MRL  function  is  decreasing  for  the  case  =  2.5.  Finally,  shape 
parameter  <  1  indicates  an  item  that  is  improving  with  time.  Note  for  the  case  =  0.5 
the  hazard  function  is  decreasing  while  the  MRL  function  is  increasing. 

MRL  and  the  Aging  Process.  Bryson  and  Siddiqui  (1969:  1472)  define  aging  as 
“the  phenomenon  whereby  an  older  system  has  a  shorter  remaining  lifetime,  in  some 
statistical  sense,  than  a  newer  or  younger  one.”  The  authors  develop  several  criteria  that 
provide  statistical  evidence  of  aging  based  on  lifetime  distribution  classes.  Among  them 
are  the  increasing  hazard  rate  (IHR)  and  DMRL  functions.  They  define  IHR  as 
h{t2 )  ^  h{t^ )  for  all  r2  ^  ^  0  and  DMRL  as  m{t2 )  <  m{t^ )  for  all  t2>t^>0 .  Recall 

plots  of  these  functions  for  the  Weibull  distribution  with  shape  parameter  P  =  2.5  are 
shown  in  Figure  1 .  “Net  decreasing  mean  residual  lifetime,”  more  commonly  called  “new 
better  than  used  in  expectation”  (NBUE),  is  another  aging  criterion  specified  by  the 
authors.  NBUE  is  related  to  the  MRL  function  and  defined  as  m(0)  >  m(t)  for  all  t>0. 
Implications  between  these  distribution  classes  are  shown  below. 

IHR  =>  DMRL  =>  NBUE 

Note  the  implications  between  the  classes  are  unidirectional.  The  DMRL  criterion  is  less 
restrictive  than  the  IHR  criterion  in  that  the  hazard  rate  may  not  be  strictly  increasing, 
while  the  MRL  function  associated  with  the  same  distribution  is  decreasing.  For  this 
reason  Chen,  Hollander,  and  Langberg  (1993:120)  conclude  “the  DMRL  class  may  be 
more  appropriate  than  the  IHR  class  when  the  underlying  physical  process  suggests  wear- 
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out  but,  the  failure  rate  is  expected  to  fluctuate  and  not  satisfy  the  IHR  criterion.” 
Similarly,  the  NBUE  criterion  is  less  restrictive  than  the  DMRL  criterion  in  that  the  MRL 
function  may  fluctuate  and  not  be  strictly  decreasing,  while  the  same  distribution  exhibits 
the  NBUE  characteristic. 

One  might  further  suppose  that  knowledge  of  the  hazard  rate  at  time  t  implies 
knowledge  of  the  MRL  at  time  t.  Specifically,  does  h(t)  =  0  imply  m(t)  =  0  and  h’(t)  <  0 
imply  m'{t)  >  0?  Muth  (1975:  15)  illustrates  the  somewhat  surprising  answer  is  no.  He 
considers  the  distribution  defined  by  the  density,  hazard,  and  MRL  functions  shown  in 
Table  3.  Plots  of  the  hazard  and  MRL  functions  for  this  distribution  are  shown  in 
Figure  2. 


Table  3.  Example  Distribution  Representations 


m 

[(l  +  2.3t")"  -4.6t]expj^-t-^t'  ^ 

hit) 

(l  +  2.3t"f -4.6t 
l  +  2.3t" 

mif) 

1 

l  +  2.3t^ 

15 

10 
hO) 

5 

0 

Figure  2.  Hazard  and  MRL  Functions  for  the  Distribution  Defined  in  Table  3 
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Note  that  hazard  function  for  this  distribution  is  not  strictly  increasing,  yet  it  does  have  a 
decreasing  MRL  function.  Furthermore,  note  that  h\t)  <  0  in  the  interval  [0,  0.35)  while 
m{t)  <  0  in  the  same  interval.  Finally,  note  that  /i'(0.35)  =  /i(0.35)  =  0  while  m(0.35)  >  0 
and  m'(0.35)  <  0.  This  example  highlights  an  important  difference  between  the  hazard 
and  MRL  functions.  Knowledge  of  the  hazard  function  at  time  t  provides  information 
only  about  the  immediate  future,  while  knowledge  of  the  MRL  function  at  time  t  provides 
information  about  the  complete  future.  Muth  highlights  this  as  another  important 
advantage  of  the  MRL  function  over  the  hazard  function  for  characterizing  the  aging 
process  of  an  item. 

It  is  well  known  that  the  exponential  distribution  possesses  the  “memoryless” 
property.  An  item  with  an  exponential  lifetime  distribution  is  always  as  good  as  new 
since,  at  any  time  t,  its  MRL  is  the  same  as  at  time  zero.  The  item  has  no  “memory”  of 
the  passage  of  time.  Muth  defines  the  concept  of  “virtual  age”  as  a  measure  of  memory. 
The  virtual  age  of  an  item  is  “the  amount  by  which  the  expected  life  has  been  reduced 
due  to  elapsed  time.”  The  virtual  age,  L.  of  an  item  can  be  expressed  in  terms  of  the 
MRL  function: 

ty{t)  =  m(0)  -  m{t).  ( 4  ) 

If  the  virtual  age  is  zero,  the  item  is  defined  has  having  “no  memory.”  Note  that  for  the 
exponential  distribution  ty  is  zero  for  all  t.  As  another  example,  consider  an  item  with  a 
Weibull  lifetime  distribution  with  shape  parameter  p  =  2.5  and  scale  parameter  a=  1.0. 
The  virtual  age  of  the  item  at  time  f  =  1 .0  is 

fv(l)=  m(0)-m(l) 

=  0.887-0.289  =  0.598. 
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This  item  does  have  memory,  but  it  is  not  “perfect”  in  that  tv(t)  <  t.  Muth  denotes  the 
special  case  of  tv(t)  =  r  as  “perfect  memory.”  Note  that  0  <  tv(t)  <  fi  for  all  t  for  an  item 
with  a  NBUE  distribution  and  tv{t2)  ^  tv(ti)  for  all  t2  >  if  an  item  is  characterized  with  a 
DMRL  function. 

Non-Parametric  Estimation  of  the  MRL  Function.  Suppose  we  desire  to 
characterize  the  aging  process  of  an  item  given  a  random  sample  from  the  item’s 
unknown  lifetime  distribution  F{t).  One  approach  to  this  problem  is  to  fit  a  parametric 
distribution  to  the  data  and  use  the  associated  hazard  and  MRL  functions  to  describe  the 
aging  process.  The  dangers  inherent  in  this  approach  are  well  known,  especially  if  the 
sample  size  is  small  (Leemis,  1995:252).  Another  approach  is  to  directly  characterize  the 
aging  process  through  non-parametric  analysis  of  the  MRL  function.  Bryson  and 
Siddiqui  (1969:  1485)  demonstrate  the  value  of  plotting  the  empirical  MRL  function  for 
this  purpose.  The  sample  MRL  function  has  advantages  over  the  sample  density  or 
hazard  functions  in  that  the  sample  MRL  function  smoothes  fluctuations  inherent  in 
observed  data  (Muth,  1975:22). 

The  sample  MRL  function,  is  formed  by  substituting  the  sample  survivor 
function,  for  S(t)  in  (2)  to  obtain 


m„(0  = 


1 

Sn(t) 


1 {u)du . 

t 


(5) 


Guess  and  Proschan  (1985:8-9)  provide  details  for  computing  the  sample  MRL  function. 
Their  results  for  the  case  of  no  ties  in  the  data  are  shown  below. 
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m„(t)  = 


i=k+l 


n  —  k 


for  and  ^  =  0, 1, ..., n-1  (6) 


where  n  =  number  of  sample  observations; 

Xi  =  i'*  ordered  observation  where  Xi.j  <  Z,  <  Z,+;; 

Zo  =  0. 


They  also  provide  details  for  computing  the  empirical  MRL  when  ties  exist  in  the  data 
not  included  here. 

Chen,  Hollander,  and  Langberg  (1983)  propose  a  sample  MRL  function  for  censored 
data  by  using  the  Kaplan-Meier  estimate  for  the  survival  function,  SKmdf)-  Let  {Zj,5j)  be 
the  ordered  set  of  observations  where  ^  =  0  if  observation  j  is  censored  and  5j  =  1  if 
observation  j  is  not  censored.  The  Kaplan  Meier  estimate  of  the  survival  function  is  then 

1  ,  0  <  t  <  Zi 


'KME 


n  i 

(o=n 


,  Zft-i  <  t  <  Z* 


(7) 


0  ,  t  >  Z„ 

Substituting  equation  (7)  for  5„(t)  in  equation  (5)  the  MRL  function  becomes 


n-1 


^KME  (0  —  ■ 


^  KME 


(0  It 


^^^KME  (Z,.,)(Z,.,-Z,)  +  S™,(Z.)(Z,-0 


i=k 


Zk-i  <t<Zk  (  8  ) 


It  is  obvious  from  equations  (6)  and  (8)  that  m„(r)  and  (t)  respectively  are  not 
smooth  and  have  nodes  at  the  sample  observations.  Kulasekera  (1991)  obtained  a 
smooth  MRL  function  from  a  smooth  estimate  of  the  distribution  function,  which  in  turn 
is  derived  from  the  kernel  estimator  of  the  density  function.  The  kernel  density  estimator 
is  constructed  similar  to  a  histogram  whereby  a  “bump”  called  the  kernel  function,  K(x), 
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is  centered  over  each  observation.  The  kernel  function  is  itself  a  density  function  that  is 
typically  chosen  to  be  symmetric  and  mound  shaped.  The  kernel  estimate  of  the  density 
function  is  formed  by  summing  the  kernel  function  “bumps.” 

In  the  case  of  complete  data,  Kulasekera  expressed  the  sample  distribution  function  in 
terms  of  the  kernel  density  estimator  to  obtain 


n  — 


t-X, 


(9) 


where 


W{t)=  jK(u)du. 


The  parameter  h  used  in  equation  ( 9 )  is  called  the  “smoothing  parameter”  or 
“bandwidth”  and  determines  the  width  of  the  kernel  function  bumps.  As  the  value  of  h  is 
increased  the  width  of  the  kernel  function  bumps  are  increased  thereby  reducing  the 
granularity  of  the  density  and  corresponding  distribution  function  estimates. 

Furthermore,  Kulasekera  observed  equation  (5)  could  be  expressed: 

t 

A- j(l-F„(M))dM 

m„(t)  = - 2 -  ,  t>0.  (10) 

1-F„(0 

Finally,  he  derived  the  kernel  estimate  of  the  MRL  function,  ms(t),  by  substituting 
equation  (9)  in  (10)  to  obtain 


X-j(l-F^(u))du 
m(t)  =  —  " 


n 


(11) 
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When  the  data  are  censored,  Kulasekera  uses  the  Kaplan  Meier  estimator  (7)  and  the 
kernel  density  estimator  to  derive  the  estimated  MRL  function  (t) .  The  derivation 
is  analogous  to  the  complete  data  case.  The  resulting  estimate  of  survival  function  is 


^KMEs  (0  —  1  ^  (^KME  ^ KME  (^y+1  )w{ 


y=i 


(12) 


Note  that  (Z^. )  -  8^,^^  )  is  merely  the  size  of  the  jump  in  the  Kaplan  Meier 

survival  function  at  observation  j.  Using  the  result  of  equation  (12)  in  equation  (10),  the 
smooth  estimator  of  the  censored  data  MRL  is 


M  t 

J  ^ KME  (w)fi?M  J"  S {u)du 


^KMEs  (0  —  ' 


'KMEs 


xt) 


(13) 


Kulasekera  finds  that  the  kernel  estimators  of  the  MRL  are  a  much  better  approximation 
over  the  empirical  MRL  function  for  certain  choices  of  the  bandwidth,  h. 

Yet  another  technique  for  computing  the  MRL  function  with  censored  data  is  derived 
from  the  Piecewise  Exponential  Estimator  (PEXE)  of  the  survivor  function  introduced  by 
Kim  and  Proschan  (1991).  The  premise  of  the  PEXE  is  to  estimate  the  average  failure 
rate  in  each  interval  between  successive  failures,  fit  an  exponential  survivor  function  to 
each  interval,  then  piece  the  interval  survivor  function  estimates  together  to  obtain  the 
system  survivor  function.  The  PEXE  estimate  of  the  hazard  function  is 

=  — i - ,  Zj<t<Zk  (14) 

X(n-0(Z,,,-Z,) 

i=j 

where  Zj  =  is  the  first  observed  failure  (^=  1)  <  t; 

Zk  =  is  the  first  observed  failure  (5*=  1)  >  t. 
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Note  equation  (14)  is  (observed  number  of  failures)  /  (observed  total  time  on  test)  for 
each  interval  between  observed  failures.  It  is  well  known  that  this  is  the  form  of  the 
maximum  likelihood  estimate  (MLE)  of  the  failure  rate  for  the  exponential  distribution. 
Let  Zj,j=\...m  be  the  ordered  set  of  observed  failure  times  where  m  <  « is  the  total 
number  of  observed  failures.  The  resulting  PEXE  survivor  function  is 


Zi-\  <  t  <  Zi. 


(15) 


An  important  observation  from  equation  (15)  is  that  there  is  no  attempt  to  extrapolate  the 
estimate  of  the  survivor  function  beyond  the  last  observed  failure  time.  The  authors  point 
out  three  key  advantages  of  the  PEXE  over  the  KME  of  S{t).  First,  SpExsit)  is  a 
continuous  function,  although  it  does  have  nodes  at  the  observed  failure  times,  whereas 
SKME{t)  is  a  step  function.  Second,  SpExsit)  is  responsive  to  the  location  of  the  censored 
observations  between  observed  failures  while  SxMEit)  is  not.  Third,  the  step  function 
nature  of  SxMEit)  tends  to  overestimate  the  underlying  survivor  function  while  SpExsit) 
does  not.  The  PEXE  estimate  of  the  MRL  function  was  proposed  by  Joe  and  Proschan 
(1981)  and  is  formed  by  substituting  SpExsit)  of  equation  (15)  in  equation  (3)  to  obtain 


^PEXE  — 


^PEXE  (0 


(16) 


Complete  Data  Tests  for  MRL  Distribution  Classes.  Suppose  non-parametric 
analysis  of  observed  data  suggests  the  underlying  distribution  is  DMRL  or  NBUE. 

There  are  several  statistical  tests  described  in  the  literature  to  formally  test  this 
observation.  The  boundary  between  DMRL  or  NBUE  (deteriorating  item)  and  an  IMRL 
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or  new  worse  than  used  in  expectation  (NWUE)  (improving  item)  is  the  exponential  life 
distribution  which  has  a  constant  MRL  function.  Therefore  we  wish  to  test: 

Hq:  the  life  distribution  is  exponential 


vs. 


Hi:  the  life  distribution  has  a  DMRL  and  is  not  exponential 


or 


H2:  the  life  distribution  is  NBUE  and  is  not  exponential. 

The  dual  classes  of  MRL  and  NWUE  are  tested  under  the  alternatives 

Hi':  the  life  distribution  has  an  MRL  and  is  not  exponential 


and 


H2':  the  life  distribution  is  NWUE  and  is  not  exponential 


respectively. 

Hollander  and  Proschan  (1975)  introduced  the  V*  test  statistic  to  test  the  DMRL 
(MRL)  alternatives.  They  developed  their  test  statistic  by  observing  that  the  quantity 

D{s,t)  =  5'(5)5(0[»x(i')-?n(0]  (17) 

is  a  weighted  measure  of  the  deviation  of  Hi  from  Ho  for  s<t.  Note  that  D{s,t)  =  0  if 
and  only  if  Ho  is  true.  An  average  measure  of  the  deviation  of  H]  from  Ho  is  found  by 


A(F)  =  Jj  Dis,  t)dF{s)dF{t). 

S<t 

Using  the  empirical  estimates  for  the  parameters  in  equation  (18),  they  derived 

^  n  i=\ 


(18) 


X 


(19) 


where 


C:„ 


4. 3  A  ‘2  /^2.  I3  I2  1.2  1. 

— I  — 4wi  +  3w  I - n  H — n - 1  H — i. 

3  2  2  2  6 


is  asymptotically  N(0, 1/210),  thus  the  DMRL  statistic  is  (210n)^'^^y*.  Large  values 
of  the  test  statistic  indicate  a  DMRL  while  small  values  indicate  an  MRL.  Although  not 
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presented  here,  the  authors  also  develop  a  K*  statistic  to  test  for  distributions  that  are 
NBUE. 

Aly  (1990)  proposed  a  test  for  NBUE  based  on  the  quantity 


7(5)  =  |5(0(l  +  ln5(t)>/t. 


(20) 


Note  that  =  0  if  Hq  is  true  and  7(>5)  ^  0  under  H2.  Equation  (20)  expressed  in  terms 
of  F  becomes 


Y{F)  =  \tj{F{t))iF{t) 


(21) 


where  /(«)  =  2  +  ln{\-u). 

The  quantity  =  y{F^  is  formed  by  substituting  the  empirical  distribution  for  F  in 
equation  (21).  This  quantity  expressed  in  computational  form  is 

1  ^ 


n 


(22) 


The  statistic  tn  IX  converges  to  a  N(0,1)  random  variable  as  the  sample  size  n 
becomes  large.  Because  the  convergence  is  slow,  a  the  author  developed  a  modified 
form  of  the  test  statistic  with  a  faster  rate  of  convergence  and  shows  it  to  be 


1"% -Kx) 

aX 


(23) 


where 


1  « 


K =i+-SH 


1- 


7-1 


n 


1  " 

-y 


|l  +  ln 

1 

«  )\ 
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Ho  is  rejected  in  favor  of  H2  for  large  values  of  r„  and  is  rejected  in  favor  of  H2'  for  small 
values  of  r„.  The  author  concludes  the  r„  statistic  outperforms  Hollander  and  Proschan’s 
V*  statistic  in  terms  of  the  Pitman  asymptotic  relative  efficiency  (ARE)  for  the  Weibull 
distribution. 

Ahmad  (1992)  proposes  a  new  test  for  DMRL  that  he  claims  is  easier  to  compute  and 
performs  better  than  Hollander  and  Proschan’s  V*  statistic.  However,  Kumazawa  (1993) 
asserts  this  test  is  equivalent  to  Hollander  and  Proschan’s  K*  statistic  for  NBUE.  The 
development  of  Ahmad’s  test  is  presented  here  since,  as  discussed  later  in  this  chapter, 
his  work  was  extended  to  the  censored  data  case. 

Ahmad  based  his  test  on  the  intuitive  notion  that  if  m(t)  is  decreasing  and 
differentiable,  then  m'{t)  <  0  for  t  >  0.  The  first  derivative  of  the  MRL  function  is  given 

by 


m'(t) 


(24) 


Note  that  m'(t)  <  0  if  the  quantity  -  /(0|  S(u)du  j  >  0.  Therefore  an  average 
measure  of  the  deviation  of  H2  from  Ho  is 


Sp  =  S(u)duyt . 


(25) 


Substituting  the  sample  survivor  function  5„  for  S,  the  symmetrized  U-statistic  form  of 
equation  (25)  becomes 


=- 


1 


n(«  — 1) 


X(3X,-Xp. 


i<j 


(26) 
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The  quantity  n^'^Un  is  asymptotically  ^(0,  1/3),  thus  the  test  statistic  is  Ho  is 

rejected  in  favor  of  H2  for  large  values  of  (3n)’^^C/„  and  is  rejected  in  favor  of  H2'  for 
small  values  of  {2>nf^Un- 

Censored  Data  Tests  for  MRL  Distribution  Classes.  The  three  tests  discussed 
thus  far  are  limited  in  that  they  only  apply  when  the  failure  data  is  complete.  Fortunately, 
there  are  several  tests  for  DMRL  and  NBUE  that  accommodate  randomly  censored  data 
described  in  the  literature.  Chen,  Hollander,  and  Langberg  (1983)  introduce  such  a  test 
for  DMRL.  They  extend  the  work  of  Hollander  and  Proschan  by  substituting  the  Kaplan- 
Meier  estimator  of  the  survival  function  (7)  in  equation  (18)  to  obtain 


A(^) = E  -in-/' 


(27) 


7=1 


7=1 


7=1 


where 


c,  =- 


'  n-7+1 

Analogous  to  the  derivation  of  the  V*  statistic,  the  Y  statistic  for  censored  data  is  then 


y"  = 


A(F) 

Ak 


(28) 


n  i-l 


where 


The  quantity  is  asymptotically  N(0,  Tq  )^  thus  the  DMRL  statistic  is  n^  ^V/To-  A 
consistent  estimator  of  V  is  shown  to  be 


^2 


n  \B,(2)  g,(3)  ^  g,(4)  ^  R,(5)  R,(6)  ^  R,(8) 

720  t;'(n-j  +  l)(n-ol  72  18  16  45  18  72^ 


=— +y  ■ 


B.(2)  S.(3)  B.(4)  8.(5)  8.(6)  8.(8) 

+  n< - 1 - 1 - ^ - 1-  ■ 


72 


18  16  45  18 


72 


(29) 
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where 


B,.(a)  =  exp 


9  i  1 ,  .  .  .  , 

When  ties  exist  in  the  data  between  censored  and  uncensored  observations,  treat  the 
uncensored  values  as  preceding  the  censored  values  when  forming  the  ordered  list  of  Z, 
values.  The  authors  caution  that  no  more  than  50%  of  the  observations  should  be 
censored  when  applying  the  test.  However,  they  exercise  the  test  in  an  example  where 
the  data  is  57%  censored  with  good  results. 

Lim  and  Koh  (1996)  extend  the  work  of  Aly  to  accommodate  randomly  censored  data 
by  using  the  Kaplan-Meier  estimate  Fkme  =  1  -  Skme  of  F  in  equation  (21)  to  obtain  the 
test  statistic 


LI  = 


7^^ khz') 


(30) 


where 


;=i 


V  A  V 

Cj  and  ft  defined  as  in  equations  (27)  and  (28)  respectively. 

The  quantity  is  asymptotically  N{0,  CT^),  thus  the  DMRL  test  statistic  is  n^'^I7Jo. 

The  authors  show  a  consistent  estimator  of  cr^  is 


Mi- 


4  ^  {n  —  i  +  l)(n  -  i) 


1  Z,.  Z 

- L.+ 


i'.2 


4  2fi  2fi 


5,(4) 


—  n 


1  Z„  Zt 
- 2.  +  . 


4  2fi  2fF 


5„(4) 


(31) 


where  Bi{a)  is  defined  as  in  equation  (29).  The  authors  find  the  Ln  statistic  compares 
favorably  with  Chen  Hollander  and  Lang’s  statistic  only  when  the  data  is  “slightly” 
censored. 
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Lim  and  Park  (1993)  offer  a  test  for  DMRL  with  randomly  censored  data  based  on  the 
work  of  Ahmad.  They  extend  Ahmad’s  work  by  substituting  the  Kaplan-Meier  estimate 
Fkme  =  1  -  Skme  of  F  in  equation  (25).  The  computational  form  of  the  resulting  test 
statistic  is 


A 1  w 


i-\ 


i~l 


.  >1  y=i 


(32) 


where  c/  is  as  defined  in  equation  (27).  The  quantity  5^  is  asymptotically  A(0,  V), 


thus  the  DMRL  test  statistic  is  n^'^5l/To.  A  consistent  estimator  of  is  shown  to  be 


n 


1  n-1 

^0  =-+S  — 

6  l)(n  -  i) 


5,.(4)-|R,(3)  +  iR,.(2) 


S,(4)-jB.(3)  +  is,(2)j 


(33) 


where  R,(a)  defined  as  in  equation  (29). 

Ho  is  rejected  for  significantly  large  (small)  values  of  n^'^Sl/Xo  in  favor  of  Hi  (H/).  At 

issue  is  whether  this  test  is  more  appropriate  for  DMRL  or  NBUE.  Although  the  authors 
claim  this  is  a  test  for  DMRL,  Ahmad’s  test  from  which  it  was  derived  was  later  proven 
to  be  a  test  for  NBUE.  Fortunately,  the  matter  is  inconsequential  for  this  research  effort 
as  both  DMRL  and  the  less  restrictive  NBUE  class  satisfy  the  aging  requirement. 


MRL  Numerical  Example 

The  purpose  of  this  section  is  to  demonstrate  the  MRL  concepts  described  above  via  a 
numerical  example.  I  used  the  Weibull  distribution  with  shape  parameter  =  1.5  and 
scale  parameter  a  =  1.0  to  generate  n  =  100  random  variates  used  as  the  “observed”  data 
set  for  the  example.  This  technique  allows  easy  comparison  of  the  non-parametric 
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analysis  to  the  true  underlying  distribution.  MathSoft’s  Mathcad  6.0  software 
(Appendix  A)  was  used  to  generate  the  random  variates,  perform  all  calculations,  and 
generate  the  plots  illustrated  in  this  section.  The  ordered  complete  data  set  is  shown  in 
Table  4.  The  example  is  considered  in  two  parts.  The  complete  data  case  is  examined 
first  followed  by  results  for  randomly  right-censored  data. 


Table  4.  Weibull  (1, 1.5)  Ordered  Data  Set 


Obs.# 

Xi 

Obs.# 

Obs.# 

Xi 

Obs.# 

Xi 

1 

0.0218 

26 

0.4442 

51 

0.7368 

76 

1.3439 

0.051 1 

27 

0.4503 

52 

0.7536 

77 

1 .3848 

0.0666 

28 

0.456 

53 

0.7576 

78 

1 .3926 

0.1049 

29 

0.46 

54 

0.7578 

79 

1 .3984 

0.1267 

30 

0.4735 

55 

0.7818 

80 

1.4511 

6 

0.1393 

31 

0.4888 

56 

0.7955 

81 

1 .4676 

0.2598 

32 

0.5033 

57 

0.8061 

82 

1 .4766 

0.2801 

33 

0.508 

58 

0.8265 

83 

1.5128 

0.3022 

34 

0.5125 

59 

0.8416 

84 

1.5209 

10 

0.3078 

35 

0.5268 

60 

0.848 

85 

1.5269 

11 

0.3091 

36 

0.5336 

61 

0.851 

86 

1.5414 

12 

0.3119 

37 

0.5601 

62 

0.8594 

87 

1 .5422 

13 

0.3129 

38 

0.6231 

63 

0.8821 

88 

1 .5634 

14 

0.3155 

39 

0.6366 

64 

0.8987 

89 

1 .5684 

15 

0.3197 

40 

0.6387 

65 

0.8993 

90 

1 .6375 

16 

0.3197 

41 

0.6404 

66 

0.901 

91 

1.6544 

17 

0.3363 

42 

0.6559 

67 

0.9856 

92 

1.6779 

18 

0.3379 

43 

0.66 

68 

1.0324 

93 

1.6955 

19 

0.3441 

44 

0.6775 

69 

1.1028 

94 

1.7888 

20 

0.3593 

45 

0.6777 

70 

1.1234 

95 

2.0166 

21 

0.3907 

46 

0.7051  , 

71 

1.1613 

96 

2.3462 

22 

0.3957 

47 

0.7103 

72 

1.1755 

97 

2.8135 

23 

0.4328 

48 

0.7191 

73 

1.1837 

98 

2.8182 

24 

0.4337 

49 

0.7251 

74 

1.2054 

99 

3.4563 

25 

0.4433 

50 

0.7363 

75 

1.2368 

100 

3.5434 

Complete  Data. 

Empirical  MRL  Functions.  Equation  (6)  was  used  to  generate  the  empirical 
MRL  function  and  a  smooth  approximation  of  the  MRL  function  was  obtained  using 
Kulasekera’s  result  (11).  The  Epanechnikov  kernel  function  (Silverman,  1986:  43) 
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\t\  <  a/5 


(34) 


3  1  ^ 

K{t)  =  —=  , 

4V5l,  5  J 

0,  otherwise 

was  used  for  the  smooth  approximation.  Figure  3  shows  a  plot  of  the  normal  and 
smoothed  empirical  MRL  functions  contrasted  with  the  true  MRL  function  of  the 
underlying  distribution.  A  smoothing  parameter  /i  =  0. 1  was  used  in  this  case.  Figure  4 
shows  the  same  plots,  but  with  smoothing  parameter  h  =  0.2  to  demonstrate  its  effect  on 
the  resulting  MRL  function  approximation. 


Smooth  MRL 
~  Sample  MRL 
■  ■  True  MRL 

Figure  3.  Normal  and  Smoothed  (h  =  0.1)  MRL  Functions 
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1  2  3 


Smooth  MRL 
~  Sample  MRL 
True  MRL 

Figure  4.  Normal  and  Smoothed  {h  =  0.2)  MRL  Functions 

The  overall  impression  of  the  empirical  MRL  functions  is  that  the  underlying 
distribution  does  not  have  a  strictly  decreasing  MRL  function,  but  may  posses  the  NBUE 
characteristic.  Note  the  empirical  MRL  functions  are  a  fair  approximation  of  the  true 

for  t  <  1.5.  When  t>  1.5,  the  estimated  MRL  functions  over  estimate  the 
true  MRL  function  in  all  cases.  Examination  of  the  sample  survivor  curve  shown  in 
Figure  5  offers  some  insight  into  this  phenomenon.  Note  that  the  sample  survivor 
function  underestimates  the  true  survivor  function  in  the  vicinity  1.5  <  t  <  2  while  the  tail 
of  the  survivor  function  is  exaggerated  for  t>2.  Therefore,  the  numerator  of  equation  (5) 
is  over  estimated  while  the  denominator  is  under  estimated  causing  the  poor  performance 
of  the  empirical  MRL  functions  in  this  area.  Also  note  that  progressive  smoothing  of  the 
empirical  MRL  function  results  in  increased  loss  of  detail,  but  can  lead  to  a  better 
approximation  if  the  fluctuations  are  suspect  not  to  be  inherent  in  the  tme  MRL  function. 
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Sample  Estimate 
True 

Figure  5.  Complete  Data  Survivor  Function 

Statistical  Tests.  Hollander  and  Proschan’s  V*  test  (16)  was  applied  to  the  data 
set  to  test  for  DMRL.  We  would  expect  the  test  to  favor  of  Hi  over  Ho  since  we  know  the 
underlying  parametric  distribution  is  DMRL.  The  resulting  test  statistic  value 
(210n)*^^y*  =  1.289  yields  a  p-value  =  0.099.  This  test  does  not  strongly  confirm  the 
DMRL  of  the  true  underlying  distribution,  but  does  agree  with  the  impression  of  the 
empirical  MRL  estimate. 

Aly’s  Tn  test  (23)  and  Ahmad’s  Un  test  (26)  were  also  applied  to  test  for  NBUE.  The 
resulting  test  statistics  were  2.926  and  4.008,  with  corresponding  p-values  of  0.0017  and 
3.07x10  ®  respectively.  These  tests  strongly  confirm  the  NBUE  characteristic  of  the 
underlying  distribution  and  the  impression  of  the  empirical  MRL  functions. 

Censored  Data.  The  data  in  Table  2  was  randomly  censored  using  the  method 
described  by  Kulasekera  (1990;  100).  Let  Ti,. . .  ,T„  be  a  random  sample  from  the 
censoring  distribution  G.  Then  the  randomly  right  censored  data  is  given  by  the  pairs 
(Z„  di)  where 

Zi  =  min  (Xi,Yi) 

di=  l  if  Xi  <  Yi  and  0  if  X,  >  Y. 
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I  used  an  exponential  censoring  distribution  with  rate  parameter  A  =  0.8  to  obtain  43% 


censoring  of  the  data  set.  Table  5  shows  the  resulting  ordered  (Z/,5,)  random  censored 
data  pairs. 


Table  5.  Censored  Data  Ordered  Pairs 


Obs.  # 

z, 

m 

Obs.# 

m 

Obs.# 

m 

Obs.# 

Zi 

1 

0.0218 

1 

26 

0.3119 

1 

51 

0.5268 

1 

76 

0.7578 

1 

2 

0.0511 

1 

27 

0.3129 

1 

52 

0.5601 

1 

77 

0.7818 

1 

3 

0.0566 

0 

28 

0.3155 

1 

53 

0.5802 

0 

78 

0.7955 

1 

4 

0.0619 

0 

29 

0.3197 

1 

54 

0.5814 

0 

79 

0.8018 

0 

5 

0.0666 

1 

30 

0.3197 

1 

55 

0.5968 

0 

80 

0.8265 

1 

6 

0.1043 

0 

31 

0.3363 

1 

56 

0.6186 

0 

81 

0.8435 

0 

7 

0.1049 

1 

32 

0.3379 

1 

57 

0.6255 

0 

82 

0.848 

1 

8 

0.1158 

0 

33 

0.3441 

1 

58 

0.6287 

0 

83 

0.851 

1 

9 

0.1267 

1 

34 

0.3593 

1 

59 

0.6387 

1 

84 

0.8561 

0 

10 

0.1383 

0 

35 

0.37 

0 

60 

0.6404 

1 

85 

0.8821 

1 

11 

0.1392 

0 

36 

0.3798 

0 

61 

0.6431 

0 

86 

0.8852 

0 

12 

0.1693 

0 

37 

0.3907 

1 

62 

0.6559 

1 

87 

0.8987 

1 

13 

0.1823 

0 

38 

0.3957 

1 

63 

0.6775 

1 

88 

0.9776 

0 

14 

0.1838 

0 

39 

0.3974 

64 

0.6777 

1 

89 

0.9937 

0 

15 

0.1961 

0 

40 

0.4318 

65 

0.6855 

0 

90 

1.0357 

0 

16 

0.2014 

0 

41 

0.4433 

1 

66 

0.6953 

0 

91 

1.1378 

0 

17 

0.2379 

0 

42 

0.4442 

1 

67 

0.7051 

1 

92 

1.1613 

1 

18 

0.2598 

1 

43 

0.4503 

1 

68 

0.7103 

1 

93 

1.1837 

1 

19 

0.2665 

0 

44 

0.456 

1 

69 

0.712 

0 

94 

1.4511 

1 

20 

0.2719 

0 

45 

0.46 

1 

70 

0.7191 

1 

95 

1.4676 

1 

21 

0.2752 

0 

46 

0.4735 

1 

71 

0.7251 

1 

96 

1.5269 

1 

22 

0.2801 

1 

47 

0.4888 

1 

72 

0.7363 

1 

97 

1.5422 

1 

23 

0.2818 

0 

48 

0.5047 

0 

73 

0.7368 

1 

98 

1.5634 

1 

24 

0.3027 

0 

49 

0.508 

1 

74 

0.7409 

0 

99 

1.7888 

1 

25 

0.3078 

1 

50 

0.5125 

1 

75 

0.7518 

0 

100 

1.8518 

0 

Empirical  MRL  Functions.  The  KME  (8),  smoothed  KME  (13),  and  the  PEXE 
(15)  of  the  MRL  function  were  applied  to  the  censored  data  set.  Plots  of  these  MRL 
function  approximations  contrasted  with  the  complete  data  MRL  function  and  the  true 
MRL  function  are  shown  in  Figure  6.  The  complete  data  MRL  function  is  included  so  an 
assessment  can  be  made  regarding  the  affect  of  the  censoring.  The  censored  data 
approximations  can  not  be  reasonably  expected  to  perform  better  than  the  complete  data 
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approximation.  A  smoothing  parameter  h  =  Q.  \  was  chosen  arbitrarily  for  the  smoothed 
KME  function  and  the  same  kernel  (24)  was  used  as  in  the  complete  data  case.  Note  that 
of  each  of  the  censored  data  empirical  MRL  functions  depart  significantly  from  the 
known  parametric  result  in  that  the  DMRL  is  greatly  exaggerated.  The  departure  from 
the  complete  data  MRL  function  is  even  more  dramatic  when  t  >  1.5. 


PEXE  MRL  Function 


Smoothed  KME  MRL  Function 


Figure  6.  Censored  Data  MRL  Estimates 


Recall  the  empirical  estimates  of  the  MRL  function  are  an  application  of  equation  (3) 
with  empirical  estimates  for  S{t).  Therefore,  examination  of  the  empirical  survival 
functions  should  offer  some  insight  to  the  decreased  performance  of  the  censored  data 
empirical  MRL  estimates.  Figure  7  shows  the  PEXE,  KME,  and  smoothed  KME 
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survivor  functions.  Note  that  each  of  these  functions  fail  to  preserve  the  “tail”  of  the  true 
survivor  function.  In  contrast  to  the  complete  data  set  where  at  least  sparse  information 
existed  in  the  tail  of  the  distribution,  the  censored  data  set  contains  no  information  for 
t>  1.9.  This  highlights  a  major  limitation  of  each  of  the  empirical  survivor  function  and 
corresponding  MRL  function  estimates  described  in  the  literature.  None  of  the  estimates 
attempt  to  extrapolate  or  estimate  the  survivor  function  past  the  last  observation.  Since 
m{t)  =  0  when  S{t)  =  0,  the  truncated  survivor  functions  cause  the  empirical  MRL 
functions  to  drop  off  faster  resulting  in  a  more  pronounced  DMRL.  This  effect  is 
exaggerated  as  either  the  amount  of  censoring  increases  or  the  sample  size  decreases.  A 
“semi-parametric”  technique  to  solve  this  problem  is  developed  in  Chapter  3. 

KME  PEXE 


Figure  7.  Censored  Data  Survivor  Function  Estimates 
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Statistical  Tests.  Chen  Hollander  and  Langberg’s  test  was  applied  to  the 
censored  data  to  check  for  DMRL.  The  resulting  test  statistic  value  =  1.086  with 

a  corresponding  p- value  =  0.139.  This  test  did  not  confirm  the  DMRL  of  the  underlying 
distribution  even  though  the  empirical  MRL  functions  suggest  a  strongly  decreasing 
MRL.  Lim  and  Koh’s  test  for  NBUE  resulted  in  a  test  statistic  =  8.407  with  a 

corresponding  p-value  «  0.001.  Lim  and  Park’s  test  for  NBUE  resulted  in  a  test  statistic 

5l/Xo  =  4.745  with  a  corresponding  p-value  <  0.001 .  These  tests  strongly  confirm 

the  NBUE  characteristic  of  the  underlying  distribution  and  the  impression  of  the 
empirical  MRL  functions. 

Preventive  Maintenance 

Policy  Descriptions.  Three  preventive  maintenance  policies  that  have  been 
extensively  covered  in  the  literature  are  the  Age  Replacement  Policy  (ARP),  the  Block 
Replacement  Policy  (ARP),  and  the  Opportunistic  Replacement  Policy  (ORP).  The  ARP 
is  common  when  considering  simple  systems  whereas  the  BRP  and  the  ORP  are 
appropriate  for  complex  systems.  The  term  “replaced”  is  synonymous  with  the 
expression  “repaired  to  as  good  as  new”  for  the  following  discussion.  Barlow  and  Hunter 
(1960)  define  the  ARP  (Policy  I)  whereby  the  system  is  replaced  upon  failure  or  after  a 
time  T  since  the  last  feiilure,  whichever  occurs  first.  In  the  limit  as  T  approaches  infinity, 
the  ARP  reduces  to  replacing  the  system  upon  failure  and  no  preventive  maintenance  is 
performed.  Figure  8  illustrates  representative  ARP  time  sequences. 

Barlow  and  Hunter  also  define  a  BRP  (Policy  II)  whereby  the  system  is  replaced  when 
the  time  t  is  a  multiple  of  the  replacement  period  T,  that  is  when  t  =  kT  where  k  = 

1,2,3,...,  regardless  of  the  number  of  intervening  failures.  Upon  failure  a  “minimal 
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repair”  is  made  whereby  only  the  failed  component  is  replaced.  After  such  a  minimal 
repair  the  authors  assume  the  system  failure  rate  is  not  significantly  disturbed  due  to  the 
aging  of  the  other  components.  In  the  limit  as  the  replacement  period  T  approaches 
infinity,  the  BRP  reduces  to  a  policy  of  replacing  individual  components  as  they  fail.  No 
preventive  maintenance  is  performed  and  the  system  is  never  returned  to  “as  good  as 
new”  condition.  Figure  9  illustrates  representative  BRP  time  sequences. 


ARP(7) 


O  Scheduled  System  Preventive  Maintenance 
X  Failure  with  complete  Repair 

Figure  8.  Age  Replacement  Policy 
BRP(7) 


O  Scheduled  System  Preventive  Maintenance 
A  Failure  with  minimal  Repair 

Figure  9.  Block  Replacement  Policy 
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The  ORP  was  introduced  by  Radner  and  Jorgenson  (1963)  and  modified  by  Gertsbakh 
(1977).  The  ORP  is  a  combination  of  the  ARP  and  BRP  policies  and  may  be  beneficial 
when  the  cost  of  replacing  the  system  is  less  than  sum  cost  of  replacing  individual 
components.  Under  the  ORP  individual  failed  components  are  replaced  if  a  failure 
occurs  within  time  t,  the  system  is  replaced  if  a  failure  occurs  within  the  interval  t  to  T, 
and  the  system  is  replaced  at  time  T  since  the  last  system  replacement  if  no  failures  have 
occurred.  Note  that  the  ORP  reduces  to  the  ARP  if  t=  0,  and  equates  to  the  BRP  if  t=  7. 
Figure  10  illustrates  representative  ORP  time  sequences. 


ORP(x,7) 


A  A — ( 

- 2^ - A-€ 

9 - A — t 

<-x-^ 

◄ - T - ► 

5 _ Q 

9 - ^ 

ORP(T,T) 

H  A 

- A — e 

■•-T-*!  j. 

◄ - T — 

r>i 

- ©-^ 

^x->(  1 

< —  T — ►( 
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'  LA 
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- 

- ±3 - ^ 

►  ◄ —  T  — ► 

ORP(0,r) 

^ - © - X - Q-  X  X - 


O  Scheduled  System  Preventive  Maintenance 
X  Failure  with  complete  Repair 
A  Failure  with  minimal  Repair 

Figure  10.  Opportunistic  Replacement  Policy 
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Of  the  three  policies,  the  BRP  is  generally  considered  the  easiest  to  implement  since 
the  preventive  maintenance  does  not  depend  on  the  system  failure  history  and  can  be 
scheduled  in  advance.  Another  advantage  of  this  policy  is  that  no  documentation  of 
failures  or  time  in  service  of  the  system  or  any  of  the  components  is  required  for  its 
implementation.  The  ARP  is  slightly  more  difficult  to  implement  than  the  BRP  in  that 
system  time  in  service  must  be  tracked  to  determine  when  preventive  maintenance  must 
be  performed  in  the  event  of  no  failures.  The  ORP  is  the  most  difficult  of  all  to 
implement  since  system  time  in  service  and  the  parameters  t  and  T  must  be  tracked  to 
determine  the  appropriate  maintenance  action  to  take.  These  considerations  are  not 
directly  considered  in  the  cost  models,  but  may  be  a  criteria  for  deciding  on  a  policy 
when  two  or  more  have  otherwise  equivalent  costs. 

Policy  Costs.  The  cost  of  implementing  each  policy  must  be  determined  in  order  for 
the  best  policy  to  be  chosen  for  a  given  situation.  Let  Cfht  the  cost  of  a  system  failure 
and  Cp  be  the  cost  of  a  planned  system  replacement  where  Cf>Cp.  The  long  run  average 
cost  per  unit  time  of  the  ARP  in  terms  of  the  replacement  period  T,  Ca{T)  is  well 
documented  in  the  literature  and  is  given  by 


c.(r)  = 


(35) 


where  SQ{t)  is  the  system  survivor  function. 

The  denominator  of  equation  (35)  is  the  expected  time  of  each  replacement  period.  This 
equation  is  easily  solved  given  5o(t).  When  the  system  lifetime  distribution  is  not 
explicitly  known,  a  non-parametric  cost  model  for  the  ARP  can  be  obtained  with  an 


33 


estimate  of  So(f)  from  failure  data.  In  the  limit  as  T  approaches  infinity  the  average  cost 
per  unit  time  of  the  ARP  becomes 

c,(“)  =  ^.  (36) 

/X 

To  determine  the  cost  of  the  BRP  and  ORP  policies,  the  cost  of  replacing  individual 
components  must  be  considered.  Spearman  (1986)  defines  the  costs  Cf  and  Cp  as  follows 

nc 

(37) 

i=\ 

nc 

(38) 

I=I 

where  c*  =  cost  (penalty)  of  a  system  breakdown 

Cs  =  setup  repair  cost 

Ci  =  cost  to  replace  component  i 

nc  =  number  of  components  in  the  system. 

Note  that  under  this  formulation  the  assumption  c/>  Cp  holds  for  Cb  >  0.  Also,  the  cost  of 
system  failure  followed  by  a  minimal  repair  of  the  failed  component  only  is  given  by 

Cfl  =  Cb  +  Cs  +  Ci.  (39) 

It  should  also  be  noted  that  the  total  cost  of  replacing  individual  components  is  greater 
than  the  cost  of  simultaneously  replacing  all  components  only  if  c,  >  0. 

Spearman  shows  the  long  run  average  cost  per  unit  time  of  the  BRP  in  terms  of  the 
replacement  period  T,  Cb(T)  is 

':,+t(';/,)M',(7') 

C,(T)  = - -  (40) 

where  W,(r)  =  E[Ni{t)]  is  the  expected  number  of  renewals  of 

component  i  by  time  t. 
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The  renewal  function  Wi{t)  must  be  determined  for  each  component  of  the  system  to 
solve  equation  (40).  W/(t)  is  found  via  the  fundamental  renewal  equation  (Barlow  and 
Proschan,  1965:  50)  shown  in  equation  (41). 

W,.  it)  =  J'  (1  +  W;.  it  -  x))dF,  ix) .  (41) 


Equation  (41)  is  difficult  to  solve  for  many  distributions,  including  the  Weibull. 
Spearman  develops  an  easy  to  implement  approximation  of  the  renewal  equation  for 
lifetime  distributions  with  unbounded  increasing  failure  rates  (including  the  Weibull 
distribution)  given  by 


where 


W(0»max{Wi(0,E(0} 

W,(0  = 


(42) 


1  +  -+  ^ 


V 


J 


When  the  component  lifetime  distributions  Fiit)  and/,(t)  are  not  explicitly  known,  they 
must  be  estimated  from  failure  data  to  apply  equation  (42).  Either  parametric  or  non- 
parametric  estimates  may  be  used.  In  the  limit  as  T  approaches  infinity  the  average  cost 
per  unit  time  of  the  BRP  becomes 


(43) 


The  long  run  average  cost  per  unit  time  of  the  ORP,  Co(t,T)  is  similar  to  the  BRP,  but 
the  denominator  of  equation  (40)  must  be  modified  since  the  average  duration  of  each 
replacement  period  is  no  longer  fixed  at  T.  Spearman  develops  the  expected  cost  of  the 
ORP  and  shows  it  to  be 
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Co(r,T)  = 


(44) 


_ 1=1 _ 

T-r 

T+  J[l-Fy  (T,M)]iM 

0 


where  Fv{to,t)  is  the  distribution  function  for  the  random  variable 

V,  the  time  until  the  next  system  renewal  given  the  system 
is  at  time  to. 

The  numerator  of  equation  (44)  is  similar  to  that  for  the  BRP  (40).  The  approach  to 
finding  a  solution  to  the  fundamental  renewal  equation  via  the  approximation  (42)  is 
unchanged.  To  solve  the  denominator,  Spearman  assumes  Weibull  component  lifetime 
distributions  and  develops  an  estimate  for  the  quantity  1-  Fvit,u): 


(45) 


where  a,  =  Weibull  scale  parameter  for  component  i; 

Pi  =  Weibull  shape  parameter  for  component  i. 

This  expression  must  be  numerically  integrated  to  calculate  the  denominator  of  the  ORP 
cost  equation. 

Maintenance  Policy  Optimization.  The  problem  of  optimizing  the  aforementioned 
maintenance  policies  has  received  extensive  attention  in  the  literature.  Almost  without 
exception,  the  literature  defines  “optimal”  in  terms  of  minimizing  cost.  The  cost 
equations  for  the  respective  policies  are  minimized  with  respect  to  the  replacement  period 
T  (t  and  T  for  the  ORP)  to  find  the  value  T*  providing  the  lowest  cost  for  each  policy. 
Once  each  policy  is  optimized,  they  can  be  directly  compared  to  determine  the  overall 
best  policy  for  a  given  situation.  This  approach  to  finding  an  optimal  preventive 
maintenance  policy  requires  the  cost  parameters  described  in  equations  (37)  and  (38)  to 
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be  well  defined.  This  is  not  always  the  case.  Consider  a  system  such  as  the  AMRAAM 
where  the  penalty  of  system  failure,  Cb,  is  degraded  mission  accomplishment  or,  in  the 
extreme  case,  loss  of  life.  In  this  case  the  cost  Cb,  while  obviously  high  warranting 
consideration  of  a  preventive  maintenance  policy,  is  subjective  and  its  relative  weight 
compared  to  the  other  cost  parameters  (expressed  in  dollars)  is  difficult  to  assess.  The 
traditional  approach  of  optimization  through  minimizing  the  cost  function  therefore  can 
not  be  directly  applied.  Instead,  we  will  attempt  to  express  system  reliability  for  each 
policy  as  a  function  of  cost.  A  direct  comparison  can  then  be  made  between  the  policies 
to  determine  which  policy  provides  the  greatest  reliability  for  a  given  cost  (or  conversely, 
given  a  specified  reliability  the  “optimal”  policy  is  the  one  with  the  lowest  associated 
cost).  The  methodology  for  this  approach  will  be  developed  in  chapter  3. 
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III.  Methodology 


Overview 

In  this  chapter  the  methodology  to  analyze  the  AMRAAM  data  set  is  defined.  The 
chapter  begins  with  a  look  at  the  MRL  function.  A  technique  to  improve  the  empirical 
MRL  functions  discussed  in  chapter  2  is  developed  and  demonstrated.  The  behaviors  of 
the  various  DMRL  /  NBUE  statistical  tests  are  also  examined  with  respect  to  various 
Weibull  shape  parameters,  sample  sizes,  and  censoring  levels.  In  the  next  section,  the 
preventive  maintenance  reliability  cost  model  is  developed  and  demonstrated.  A  cost 
parameter  sensitivity  analysis  is  performed  with  the  model  to  determine  the  best  policy 
with  respect  to  parameter  ratios.  Finally  the  performance  of  the  empirical  reliability  cost 
model  is  assessed  with  a  numerical  example. 

Mean  Residual  Life 

Semi-parametric  MRL  Function.  The  empirieal  MRL  function  estimates 
discussed  in  Chapter  2  are  an  application  of  equation  (3)  with  empirieal  estimates  for  S(t). 
However,  the  empirical  survivor  functions  generally  fail  to  preserve  the  “tail”  of  the 
distribution  as  the  amount  of  eensoring  increases.  The  truneated  survivor  functions  result 
in  empirical  MRL  function  estimates  with  an  exaggerated  DMRL.  This  effect  is 
unacceptable  for  two  reasons.  First,  we  are  attempting  to  characterize  the  aging  process 
to  determine  if  a  preventive  maintenance  policy  is  applicable.  An  exaggerated  MRL 
function  may  falsely  indicate  that  a  preventive  maintenance  policy  is  warranted  when,  in 
fact,  the  MRL  of  the  underlying  distribution  may  not  support  this  conelusion.  Second,  an 
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accurate  MRL  function  is  an  integral  piece  of  the  reliability  cost  model  developed  later  in 
this  chapter. 

A  method  is  required  to  extrapolate  the  empirical  survivor  functions,  past  the 
last  observation  in  order  to  improve  the  empirical  MRL  estimates.  We  propose  a  “semi- 
parametric”  technique  whereby  a  distribution  is  fitted  to  the  data  and  the  corresponding 

parametric  survivor  function,  S(t) ,  is  used  to  estimate  the  survivor  function  beyond  the 

/V 

last  observation.  The  general  form  of  semi-parametric  survivor  function,  (t)  is 

5„(0  =  S„(t)  ,  t<Z„  (46) 

~  S (t)  ,  t>  Zn- 


The  general  form  of  the  semi-parametric  MRL  function  m„  (t)  is  then 


(0 


1 

Snit) 


J  (u)du  +  C 

\  ‘  ) 


t<Zn 


(47) 


where  C  =  J .S {u)du . 

z„ 

The  constant  C  in  the  semi-parametric  MRL  function  accounts  for  the  area  in  the  tail  of 
the  empirical  survivor  function  that  may  be  lost  due  to  censoring.  Note  that  no  attempt  is 
made  to  estimate  the  MRL  function  beyond  the  last  observation.  If  an  estimate  is 
required  in  this  region,  the  parametric  MRL  function  corresponding  to  the  fitted 
distribution  is  used.  This  technique  is  extremely  flexible  in  that  the  constant  C  is  a 
function  of  the  last  observation,  Z„.  If  Z„  is  sufficiently  large  such  that  the  tail  of  the 
empirical  survivor  function  is  preserved  (or  perhaps  exaggerated),  C  will  be  negligible 
and  the  empirieal  MRL  function  is  not  significantly  altered.  If  Z„  is  such  that  significant 
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area  in  the  tail  of  the  survivor  function  has  been  lost,  C  increases  accordingly  and  the 
empirical  MRL  function  is  improved. 

We  demonstrate  the  semi-parametric  MRL  function  technique  with  an  extension  of  the 
numerical  example  presented  in  Chapter  2.  Weibull  distributions  were  fitted  to  the 
complete  data  (Table  4)  and  censored  data  (Table  5)  to  compute  the  constant  C  of 
equation  (47)  for  each  case.  Leemis  (1995)  provides  details  for  computing  Weibull  MLE 
shape  (p)  and  scale  (a)  parameters  applied  here.  The  shape  parameter  is  found  via  a 
Newton-Raphson  iterative  procedure: 


Pm  ~  Pi 


8(Pi) 

8\Pi) 


(48) 


r^Zf  InZ, 

where  g  (^)  =  -^  +  ^  In  Z,. - 

P  ieU  V  yP 


i=l 


8'(P)  =  - 


P 


^  r  n  V 


Sz/ 


i=i  X  '■=’  y  V  '■=' 


Y 

y 


V  '■=>  J 

r  =  number  of  observed  failures; 

U  =  set  of  observed  failures. 

An  initial  guess  Po=  1-5  was  used  to  start  the  procedure  and  a  tolerance 
£  =  0.01  >  \Pi+i-Pi\  was  used  for  the  final  result.  Once  p  is  found,  the  scale  parameter  is 
found  via 


/ 


a  = 


Izf 

V  '■=’  y 


(49) 
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The  complete  and  censored  data  cases  of  this  example  are  considered  separately. 

Complete  Data.  Recall  that  the  tail  of  the  non-parametric  survivor  function  was 
exaggerated  (Figure  5)  causing  a  significant  departure  of  the  non-parametric  MRL 
functions  from  the  true  distribution  (Figure  3).  The  semi-parametric  MRL  technique 
discussed  above  should  show  little  difference  with  the  non-parametric  results  since  the 
tail  of  the  survivor  function  was  preserved  in  the  non-parametric  result.  Application  of 
equation  (48)  resulted  in  Weibull  parameter  estimates  of  =  1.428  and  a„  =  0.993.  The 
parametric  estimates  are  reasonable  compared  to  the  true  parameters  =1.5  and  a  =  1 .0. 
The  complete  data  non-parametric  and  parametric  survivor  functions  are  shown  in 
Figure  11. 


Non-parametric 

Parametric 

Figure  11.  Complete  Data  Survivor  Function  Estimates 

Note  the  tail  of  the  parametric  estimate  is  very  small  past  the  last  observation.  In  fact,  the 
resulting  constant  C  =  0.000477.  As  desired,  the  semi-parametric  normal  and  smoothed 
MRL  functions  shown  in  Figure  12  are  indistinguishable  from  the  non-parametric 
functions  shown  in  Figure  3. 
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Smooth  MRL 
”  Sample  MRL 
'  “  True  MRL 

Figure  12.  Complete  Data  Semi-parametric  MRL  Functions 

Censored  Data.  The  three  censored  data  MRL  estimation  techniques  greatly 
exaggerated  the  DMRL  characteristic  of  the  underlying  distribution  (Figure  6)  due  to  the 
truncated  nature  of  the  associated  non-parametric  survivor  functions.  The  semi- 
parametric  MRL  technique  should  alleviate  this  problem.  Application  of  equation  (48)  to 
the  censored  data  set  resulted  in  Weibull  parameter  estimates  of  =  1.724  and  On  = 

1.082  and  the  constant  from  equation  (47)  C  =  0.01 1.  The  censored  data  KME  survivor 
function,  parametric  survivor  function,  and  the  true  survivor  function  are  shown  in  Figure 
13.  The  parametric  survivor  function  increases  area  in  the  tail  of  the  distribution  over 
that  of  the  non-parametric  result,  but  underestimates  the  area  in  the  tail  of  the  true 
survivor  function.  Still,  the  semi-parametric  MRL  functions  shown  in  Figure  14  are  a 
marked  improvement  over  the  non-parametric  functions  (Figure  6).  Note  the  semi- 
parametric  results  underestimate  the  true  MRL,  but  the  slope  of  the  DMRL  is  generally 
preserved. 
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Figure  13.  Non-parametric  and  Parametric  Survivor  Functions 


PEXE  MRL  Function 


KME  MRL  Function 


Smoothed  KME  MRL  Function 


Censored  Data  Estimate 
“  “  Complete  Data  Estimate 
True  MRL 


Figure  14.  Semi-parametric  MRL  Functions 
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MRL  Test  Performance.  The  behaviors  of  the  power  and  level  of  the  tests 
described  in  Chapter  2  with  respect  to  sample  size,  amount  of  censoring,  and  MRL 
characteristic  of  the  underlying  distribution  must  be  understood  before  accepting  the 
implied  conclusions  of  the  tests.  The  tests  described  in  Chapter  2  are  examined  with 
respect  to  these  parameters  in  this  section.  The  complete  data  and  censored  data  tests  are 
considered  separately. 

Complete  Data.  The  behaviors  of  the  V*,  Tn,  and  Un  test  statistics  with  changes 
in  the  Weibull  shape  parameter  of  the  underlying  true  distribution  were  assessed  via  a 
Monte  Carlo  simulation.  (See  Appendix  B  for  simulation  Mathcad  code.)  Weibull 
distribution  “data  sets”  were  generated  for  P  =  0.9  to  1.5  in  0.1  increments  with  thirty 
replications  performed  at  each  value  of  j8.  P- values  for  each  of  the  three  tests  were  then 
computed  at  each  replication  and  jS  value.  The  resulting  p-value  verses  P  scatter  plots  are 
shown  in  Figure  15.  The  p-value  mean  at  each  P  is  indicated  in  the  plots  to  give  a  sense 
of  the  overall  trend  of  the  tests  as  the  shape  parameter  is  varied.  Plots  for  n  =  100  and  n  = 
200  are  shown  side  by  side  to  capture  the  effect  of  sample  size  on  test  performance. 

Figure  15  suggests  the  variability  of  the  p-value  responses  is  less  for  Aly’s  T„  and 
Ahmad’s  t/„  tests  than  for  Hollander  and  Proschan’s  V*  test.  This  confirms  the  claims  by 
Aly  and  Ahmad  that  their  tests  are  more  efficient  than  the  V*  test.  All  three  tests  have  a 
large  amount  of  variability  when  the  underlying  distribution  is  exponential  (j8  =  1). 
However,  Ho  is  incorrectly  rejected  at  the  0.1  level  for  no  more  than  10%  of  the 
replications  in  this  case.  The  figure  also  suggests  the  intuitive  conclusion  that  the 
performances  of  the  three  tests  improve  as  the  sample  size  increases.  Note  the  variability 
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of  each  test  is  reduced  when  the  sample  size  is  increased  to  200  for  all  values  of  p  except 

P=  1.0. 


n  =  too  n  =  200 


mean 


Figure  15.  Complete  Data  Test  Performance 
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Censored  Data.  The  same  Monte  Carlo  simulation  procedure  used  to  analyze  the 
complete  data  tests  was  applied  to  the  censored  data  V^,  L/,  and  5n  tests  described  in 
Chapter  2.  The  resulting  p-value  versus  ^  scatter  plots  for  45%  data  censoring  are  shown 
in  Figure  16.  The  p-value  mean  is  included  in  the  plots  to  give  a  sense  of  the  overall 
trend  of  the  tests  as  the  shape  parameter  is  varied.  Plots  for  n  =  100  and  n  =  200  are 
shown  side  by  side  to  capture  the  effect  of  sample  size  on  test  performance.  The  scatter 
plots  for  Lim  and  Koh’s  L/  test  are  not  included  since  p-values  were  «  0.001  for  all  )S 
values  and  sample  sizes. 


n  =  100 


n  =  200 


mean 


Figure  16.  Censored  Data  Test  Performance  with  Sample  Size  (45%  Censoring) 


46 


The  overall  impression  from  the  scatter  plots  in  Figure  16  is  that,  as  with  the  complete 
data  tests,  the  performances  of  the  censored  data  tests  improve  with  increasing  sample 
size.  Furthermore,  the  figure  confirms  Lim  and  Park’s  assertion  that  their  5n  test  is  more 
efficient  than  Chen  Hollander  and  Langberg’s  Y  test.  Also  note  the  test  is  unreliable 
when  the  data  is  subject  to  45%  censoring  with  a  sample  size  of  n  =  100. 

Figure  17  shows  p-value  versus  p  scatter  plots  for  n  =  100  with  plots  for  9%  and  45% 
censoring  shown  side  by  side  to  capture  the  effect  of  censoring  percentage  on  test 
performance.  Again  the  scatter  plots  for  Lim  and  Koh’s  L/  test  are  not  included  since  p- 
values  were  «  0.001  for  all  cases.  The  tests  improve  dramatically  with  decreased 
censoring.  Note  the  variability  in  the  responses  decreases,  the  probability  of  rejecting  Ho 
when  Ho  is  false  increases,  and  the  probability  of  rejecting  Ho  when  Ho  is  true  decreases 
with  the  lower  censoring  percentage. 

The  extremely  low  p- values  obtained  in  all  cases  with  Lim  and  Koh’s  Ln  are  notable. 
This  phenomenon  is  consistent  with  the  results  the  authors  obtained  when  they  applied 
their  test  to  failure  data  consisting  of  21 1  observations  subjected  to  57%  random  right 
censoring.  In  this  example  they  obtained  a  test  statistic  value  =8.75  with  a 

corresponding  p-value  «  0.001.  Conversely,  Chen  Hollander  and  Langberg  obtained  a 
test  statistic  value  n^'^Y/Zo  =  L52  with  a  corresponding  p-value  =  0.064  and  Lim  and 
Park  obtained  test  statistic  value  n^'^5l/Zo  =  3.870  with  a  corresponding  p-value  <  0.001 

when  their  tests  were  applied  to  the  same  data  set.  In  any  case,  Lim  and  Koh’s  Ln  test 
was  deemed  unreliable  for  this  research  effort  since  it  never  correctly  failed  to  reject  Ho 
under  the  conditions  of  the  simulation  efforts  just  presented. 
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In  general,  the  performance  of  Chen  Hollander  and  Langberg’s  and  Lim  and  Park’s 
tests  improve  as  either  the  sample  size  increases  and/or  the  amount  of  censoring 
decreases.  Since  it  is  impractical  to  attempt  quantify  the  goodness  of  each  test  for  every 
sample  size  and  censoring  percentage  combination,  the  tests  should  be  subject  to 
simulations  similar  to  the  ones  presented  here  before  accepting  the  conclusions  of  the 
tests.  Recall  the  current  AMRAAM  data  set  consists  of  815  observations  subject  to  74% 
censoring  with  future  data  expected  to  consist  of  over  2500  observations  subject  to  80% 
censoring.  The  and  5n  were  subjected  to  the  same  Monte  Carlo  simulation  procedure 
described  above  to  assess  their  suitability  under  these  conditions  with  the  results  shown 
in  Figure  18  and  Figure  19  respectively.  The  figures  suggest  Chen,  Hollander,  and 
Langberg’s  test  is  unreliable  with  these  combinations  of  sample  size  and  censoring 
percentage.  Note  the  p- value  ~  0. 1  for  all  values  with  no  decreasing  p-value  trend  as  /? 
increases.  Lim  and  Park’s  8n  test  performs  better  in  that  there  is  a  decreasing  p-value 
trend  with  increasing  j8.  However,  the  figure  suggests  a  high  probability  of  rejecting  Ho 
when  Ho  is  true  ()8  <  1.0). 
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Figure  18.  V  and  8n  Test  Performance  with  n  =  815  and  15.5%  Censoring 
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mean 


Figure  19.  and  5n  Test  Performance  with  n  =  2500  and  80.7%  Censoring 
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Preventive  Maintenance 


Reliability  Cost  Model  Development  As  stated  in  Chapter  2,  the  traditional 
approach  of  maintenance  policy  optimization  through  minimizing  the  cost  function  can 
not  be  directly  applied  when  the  penalty  of  system  failure  is  subjective.  In  this  case,  a 
cost  model  to  express  system  reliability  as  a  function  of  cost  for  each  preventive 
maintenance  policy  is  more  appropriate.  The  cornerstone  of  this  “reliability  cost  model” 
is  the  ability  to  compute  the  costs  of  the  respective  preventive  maintenance  policies  given 
a  reliability  goal.  Recall  that  the  age,  block,  and  opportunistic  preventive  maintenance 
policies  described  in  Chapter  2  have  planned  system  renewals  at  the  replacement  period 
T.  Thus,  T  is  directly  related  to  system  reliability.  A  small  value  of  T  results  in  greater 
reliability  since  the  system  is  preventively  returned  to  new  condition  more  often,  thereby 
lowering  the  probability  of  failure.  Since  the  respective  policy  costs  are  a  function  of  T, 
a  link  between  system  reliability  and  policy  costs  is  established. 

Reliability  Determination.  The  survivor  function  of  an  item,  S{t)  =1-  F{f),  is 
the  probability  the  lifetime  of  an  item  is  greater  than  t  and  therefore  is  a  measure  of  the 
item’s  reliability  at  time  t.  One  might  suppose  that  selecting  T  such  that  S(T)  equals  the 
desired  reliability  goal  is  an  appropriate  method  for  specifying  system  reliability  with  a 
preventive  maintenance  policy.  However,  this  intuition  is  erroneous.  Consider  the 
example  of  an  item  with  an  exponential  lifetime  distribution  with  associated  survivor 
function  as  shown  in  Figure  20.  It  appears  a  preventive  maintenance  policy  can  be 
applied  with  a  value  of  T  selected  to  meet  a  reliability  goal.  We  know,  however,  a 
preventive  maintenance  policy  can  not  improve  reliability  since  the  item  is  always  “as 
good  as  new”  due  to  the  memoryless  property  of  the  exponential  distribution. 
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Figure  20.  Exponential  Distribution  Survivor  Function 


The  conflict  arises  since  the  survivor  function  is  a  measure  of  projected  reliability  starting 
from  t  =  0,  and  does  not  consider  the  probability  of  survival  at  time  t  given  the  item  has 
survived  to  some  time  between  0  and  t.  The  “virtual  age”  of  an  item  described  in  Chapter 
2  is  a  measure  of  this  conditional  passage  of  time.  Let  Ty  denote  the  “virtual  replacement 
period”  determined  from  the  survivor  function  as  in  Figure  20.  From  equation  (4),  Ty  can 
be  expressed 

Ty  =  miO)-m(T) .  (50) 


Rearranging  we  obtain 


miT)  =  m(0)-Ty.  (51) 

Given  a  value  of  Ty  determined  from  the  survivor  function,  the  actual  replacement  period 
T  can  then  be  determined  from  the  MRL  function.  This  process  applied  to  the 
exponential  item  in  the  example  is  shown  graphically  in  Figure  21.  Note  there  is  no 
solution  to  equation  (51)  given  a  value  of  7^,  >  0  since  m{t)  =  m(0)  is  constant  for  all  t. 
This  confirms  the  notion  that  a  preventive  maintenance  policy  is  not  suitable  for  an  item 
with  an  exponential  lifetime  distribution. 
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Exponential  Survivor  Function  Exponential  MRL  Function 


To  illustrate  the  process  with  an  aging  item,  consider  a  Weibull  lifetime  distribution 
with  survivor  and  MRL  functions  as  shown  in  Figure  22.  A  reliability  goal  of  0.8  yields 
a  virtual  replacement  period  =  0.5.  Applying  equation  (51)  we  obtain  m{T)  =  m(0)  -  Ty 
=  0.9  -  0.5  =  0.4.  From  the  MRL  function  the  true  replacement  period  T  is  determined  to 
be  0.7. 

Weibull(  1,2.5)  Survivor  Function  Weibull(  1,2.5)  MRL  Function 


m(7)  =  m(0)  -Ty  =  0.9  -  0.5  =  0.4 

Figure  22.  Replacement  Period  determination  with  the  Weibull  (1, 2.5)  Distribution 
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Since  it  is  somewhat  awkward  to  enter  the  survivor  function  and  MRL  function  graphs 
with  the  respective  dependant  variables,  a  slightly  modified  procedure  to  apply  the 
reliability  cost  model  is  summarized  as  follows: 

1)  Specify  a  replacement  period  T 

2)  Compute  preventive  maintenance  policy  costs  given  T 

3)  Determine  the  “virtual  replacement  period”.  Tv  =  m(0)  -  m(T) 

4)  Determine  the  system  reliability  =  S(Tv). 

This  procedure  reverses  the  direction  of  the  arrows  shown  in  Figure  22. 

Policy  Costs.  The  policy  costs  discussed  in  Chapter  2  must  be  modified  for  the 
reliability  cost  model.  Recall  the  policy  costs  for  the  age,  block,  and  opportunistic 
preventive  maintenance  policies  (equations  (35),  (40),  and  (44)  respectively)  included  the 
cost  of  a  system  failure  cj  as  defined  in  equation  (37).  However,  the  penalty  of  a  system 
breakdown  (cb)  is  not  explicitly  included  in  the  reliability  cost  model.  Therefore  the  cost 
of  a  system  failure  becomes 

nc 

Cf=Cs+'^Ci.  (52) 

i=l 

Note  that  Cf  is  now  equivalent  to  the  cost  of  a  planned  system  replacement  (cp)  as  defined 
in  equation  (38).  The  reliability  model  policy  costs  for  the  ARP,  BRP,  and  ORP 
normalized  with  respect  to  Cp  are  then  given  by 

=  ^ - ,  (53) 
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Csr(T)  = 


v 

X(c,  +  c>,(r) 

ILl _ _ _ 
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(54) 
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1=1 
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T+  J[l-F^,(T,M)]^iM 


(55) 
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Reliability  Cost  Model  Example.  In  this  section  the  complete  reliability  cost  model 
is  demonstrated  with  a  numerical  example.  A  system  composed  of  four  components  is 
considered.  The  components  are  assumed  to  have  Weibull  lifetime  distributions  with 
parameters  as  specified  in  Table  6.  A  setup  repair  cost  =  0.1  was  chosen  arbitrarily  and 
component  replacement  costs  were  specified  as  shown  in  Table  7.  From  equation  (38) 
the  cost  of  a  complete  system  replacement  is  then  Cp  =■  1.0. 


Table  6.  Component  Weibull  Distribution  Parameters 


Component 

a 

P 

1 

0.0003034 

1.2 

2 

0.0002716 

1.3 

3 

0.0002848 

1.4 

4 

0.0002736 

1.5 

Table  7.  Component  Replacement  Costs 


Component 

Ci 

^fi  “  "I" 

1 

.3 

.4 

2 

.25 

.35 

3 

.2 

.3 

4 

.15 

.25 

The  system  survivor  function  So{t)  from  equation  (1)  is 

( 

50(0  =  exp 


C  *  »  ^ 

-IkO" 


V  <•=! 


and  the  resulting  system  MRL  function  mo{t)  is  given  by 


m. 


1  ” 

(t)  =  ——jSo(u)du. 


(56) 


(57) 
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Plots  of  the  system  survivor  and  MRL  functions  are  shown  in  Figure  23. 

System  Survivor  Function  System  MRL  Function 


Figure  23.  System  Survivor  and  MRL  Functions 

The  reliability  cost  model  summarized  above  was  implemented  with  the  replacement 
period  T  varied  from  300  to  4500  hours  in  100  hour  increments.  Equation  (42)  was  used 
to  solve  the  expression  in  the  numerator  of  the  BRP  (54)  and  ORP  (55)  cost  equations. 
Equation  (45)  was  used  to  solve  the  expression  in  the  denominator  of  the  ORP  cost 
equation  (55).  An  iterative  procedure  whereby  the  parameter  t  was  varied  from  0  to  T  in 
100  hour  increments  was  used  to  optimize  the  ORP  (55)  with  respect  to  rat  each  value  of 
T.  All  computations  and  plots  were  computed  via  Mathsoft’s  Mathcad  6.0  software  on  a 
PC  (Appendix  C).  Figure  24  shows  the  resulting  policy  costs  as  a  function  of  the 
computed  reliability  at  each  value  of  T.  Figure  25  shows  the  parameters  T  and  t  plotted 
against  the  computed  reliability. 

Some  important  observations  can  be  made  from  this  example.  First,  the  policy  costs 
are  increasing  functions  with  respect  to  the  reliability  goal,  and  the  increases  are  dramatic 
when  the  specified  reliability  goal  exceeds  0.85.  Second,  the  most  cost  efficient  policy 
(in  this  case  the  ORP)  remains  the  same  regardless  of  the  specified  reliability  goal. 

Third,  the  distinction  between  the  costs  of  the  different  policies  diminishes  as  the 
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specified  reliability  goal  increases.  These  observations  are  consistent  with  Spearman’s 
results  for  the  traditional  cost  optimization  model  employed  with  several  different  cost 
coefficients  and  Weibull  distribution  parameters. 


-  BRP 
ORP 


Figure  24.  Policy  Costs  vs.  Reliability  Goal 


"  “  Tau 

Figure  25.  T and  tvs.  Reliability  Goal 


Cost  Parameter  Sensitivity  Analysis. 
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Spearman’s  results  differed  from  those  presented  in  the  example  above  in  that  he  found 
the  ARP  to  be  most  cost  efficient  policy.  He  also  found  the  ORP  to  be  degenerate  (t  =  0) 
with  the  ARP.  The  disparities  are  accounted  for  by  the  choice  of  the  various  cost 
coefficients,  notably  the  value  of  the  setup  repair  cost  (c^)  relative  to  the  sum  cost  of  the 
individual  components  (Zed-  If  Cs  is  very  large,  as  in  Spearman’s  examples,  the  ARP  is 
attractive  since  the  cost  Cs  is  incurred  less  often  than  with  the  BRP.  If  c.,  is  relatively 
small,  as  in  the  example  presented  above,  the  BRP  is  attractive  since  good  components 
are  replaced  less  often  than  with  the  ARP. 

At  issue  is  the  value  of  Cs  that  represents  the  boundary  between  the  ARP  and  the  BRP 
costs.  Figure  26  shows  a  plot  of  the  three  policy  costs  verses  the  cost  ratio  cJZci.  The 
parameter  Cs  was  varied  from  0  to  4.5  while  the  other  costs  were  held  constant  as 
described  in  Table  7  to  obtain  a  range  of  the  cost  ratio  cJZci  from  0  to  5.0.  Note  the 
boundary  between  the  ARP  and  the  BRP  is  cJZci  =  0.75.  Furthermore,  note  that  the 
ORP  is  always  at  least  as  cost  efficient  as  either  the  ARP  or  the  BRP.  As  the  ratio  cJZci 
becomes  very  large  the  optimal  value  of  T  decreases  and  the  ORP  becomes  degenerate 
with  the  ARP  as  shown  in  Figure  27.  As  the  ratio  cjZci  decreases  to  zero  the  optimal 
value  of  T  approaches  T  and  the  ORP  becomes  degenerate  with  the  BRP.  The  ORP  offers 
the  greatest  advantage  in  cost  efficiency  when  the  cost  ratio  cJZci  =  0.75. 


59 


BRP 

ORP 


Figure  26.  Policy  Costs  vs.  q/Eq  Cost  Ratio,  T=  1800 


Cost  Ratio 


Figure  27.  ORP  tvs.  cJYjCi  cost  ratio 

Empirical  Reliability  Cost  Model.  The  performance  of  the  reliability  cost  model 
applied  to  empirical  data  was  assessed  via  a  Monte  Carlo  simulation  (Appendix  C). 
Weibull  distribution  “data  sets”  were  generated  for  a  four  component  series  system  in 
accordance  with  the  parameters  in  Table  6.  The  data  set  was  censored  in  accordance  with 
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the  following  procedure.  Let  Xy  represent  the  failure  time  for  the  /*  component,  then 

the  censored  system  level  failure  data  is 

Z,  =  min  {Xa,  Xa,  Xb,  Zm,  F,  }  (58) 

where  F,  is  an  observation  from  the  censoring  distribution  G. 

I  used  an  exponential  censoring  distribution  with  rate  A  =  0.0018  to  obtain  an  average  of 
74%  random  right  censoring  of  the  system  level  data  set.  Note  that  if  a  system  level 
observation  is  censored  (Z/  =  F,),  then  the  observations  for  each  of  the  components  are 
also  censored.  If  a  system  level  observation  is  an  observed  failure  (Z,  >  F,),  then  all  but 
the  minimum  component  failure  time  are  censored.  Therefore,  component  level  data  is 
always  subject  to  greater  censoring  than  system  level  data  with  this  censoring  scheme. 

Reliability  Determination.  The  system  reliability  S(Tv)  is  a  random  variable  in 
the  empirical  model  since  it  is  a  function  of  values  calculated  empirical  estimates  of  the 
survivor  and  MRL  functions.  To  assess  the  performance  of  this  random  variable,  30 
replications  of  n  =  815  data  points  were  generated  as  stated  above.  The  system  reliability 
was  computed  for  values  of  T  from  500  to  2000  hours  in  300  hour  increments  for  each 
replication.  This  procedure  was  repeated  for  each  survivor  and  MRL  function  estimation 
technique:  KME,  PEXE,  and  smoothed  KME  respectively. 

Figure  28  shows  the  resulting  reliability  verses  T  scatter  plot  using  the  semi- 
parametric  KME  estimate  for  the  MRL  function.  The  mean  and  its  associated  95% 
confidence  interval  using  the  student’s  t  statistic  are  also  shown  to  provide  a  sense  of  the 
overall  trend  and  variability  of  the  responses.  The  true  reliability  curve  is  shown  for 
comparison.  The  95%  confidence  interval  covers  the  true  result  for  all  values  of  T 
considered  (except  T=  1700)  suggesting  no  statistical  difference  between  the  empirical 
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mean  and  the  true  result.  The  same  simulation  was  performed  using  the  non-parametric 
KME  MRL  function  with  the  results  shown  in  Figure  29.  The  random  number  seeds  used 
for  the  previous  simulation  were  repeated,  therefore  the  differences  between  the  two 
results  are  attributed  solely  to  the  different  MRL  estimation  techniques.  Note  that  the 
reliability  for  a  given  value  of  T  is  generally  underestimated  when  the  non-parametric 
MRL  function  is  used.  The  95%  confidence  interval  does  not  cover  the  true  result 
indicating  statistical  significance  between  the  empirical  and  true  result.  This 
characteristic  is  attributed  to  the  generally  exaggerated  DMRL  functions  of  this 
technique. 


Data 


Mean 

Lower  95%  Cl 
-o-  Upper  95%  Cl 
—  True 

Figure  28.  Semi-parametric  KME  Reliability  vs.  T 
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Data 

Mean 

-<»-  Lower  95%  Cl 
-o-  Upper  95%  Cl 
“  True 

Figure  29.  Non-parametric  KME  Reliability  vs.  T 


Although  the  semi-parametric  KME  result  performed  well  on  average,  the  variability 
in  the  responses  is  of  concern,  especially  for  large  values  of  T.  Figure  30  and  Figure  3 1 
show  similar  reliability  verses  T  scatter  plots  with  using  the  semi-parametric  PEXE  and 
smoothed  KME  {h  =  100)  MRL  functions  respectfully.  The  random  number  seeds  used 
for  the  KME  simulations  were  repeated  allowing  direct  comparison  of  the  results.  The 
PEXE  MRL  function  results  are  similar  to  those  for  the  semi-parametric  KME.  Although 
a  slight  improvement  in  the  mean  response  for  the  larger  values  of  T  is  detectable,  there  is 
no  significant  reduction  in  the  variability  of  the  responses. 
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Data 

Mean 

Lower  95%  Cl 
-o-  Upper  95%  Cl 
~  True 

Figure  30.  Semi-parametric  PEXE  Reliability  vs.  T 

The  smoothed  KME  MRL  function  (h  =  100)  provided  a  slight  reduction  in  the 
variability  of  the  responses  over  those  of  the  semi-parametiic  KME  and  PEXE  MRL 
functions.  A  smoothing  parameter  h  =  200  was  attempted  with  the  results  shown  in 
Figure  32.  Although  the  variability  of  the  responses  is  reduced  over  the  previous  cases, 
the  reliability  is  underestimated  for  small  values  of  T. 

In  summary,  the  semi-parametric  smoothed  KME  MRL  function  provides  the  best 
results  for  reliability  determination  of  the  three  MRL  estimation  techniques.  However, 
the  smoothing  parameter  h  must  be  chosen  judiciously.  The  smoothing  parameter  should 
be  chosen  as  large  as  possible  without  so  distorting  the  empirical  survivor  and  MRL 
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functions  that  an  accurate  reliability  determination  is  not  possible.  Also,  regardless  of  the 
empirical  survivor  /  MRL  function  estimation  technique,  the  Monte  Carlo  simulations 
show  an  accurate  estimation  of  reliability  can  not  be  assumed  when  T  is  large  and 
approaches  the  value  of  the  last  observation.  This  is  not  considered  to  be  a  significant 
limitation,  however,  since  in  most  cases  where  a  preventive  maintenance  policy  is 
considered  small  values  of  T  are  desired. 


’<xx  Data 
Mean 

Lower  95%  Cl 
-o*  Upper  95%  Cl 
True 

Figure  31.  Semi-parametric  Smoothed  ICME  Qi  =  100)  Reliability  vs.  T 
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Data 
■®“  Mean 
-<>■  Lower  95%  Cl 
-®-  Upper  95%  Cl 
“  True 

Figure  32.  Semi-parametric  Smoothed  KME  (h  =  200)  Reliability  vs.  T 


Policy  Costs.  The  costs  of  the  three  preventive  maintenance  policies  are  also 
random  variables  in  the  empirical  model  since  their  calculation  requires  estimates  of 
either  the  system  (ARP)  or  component  (BRP  and  ORP)  lifetime  distribution  functions. 
The  performances  of  the  empirical  policy  costs  were  assessed  with  the  same  Monte  Carlo 
simulation  procedure  described  above  for  the  reliability  determination.  The  empirical 
ARP  cost  was  determined  via  equation  (53)  with  KME  of  the  system  survivor  function 
used  in  place  of  Soit).  The  empirical  BRP  cost  was  determined  via  equation  (54).  The 
component  renewal  functions,  l^(t),  were  solved  via  equation  (42).  Weibull  MLE 
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parameters  were  determined  via  equation  (48)  and  used  to  compute  parametric  estimates 
for  jj,,  and  F(t)  using  the  equations  presented  in  Table  2.  The  empirical  ORP  cost 

was  determined  via  equation  (55).  The  component  renewal  functions  were  computed  as 
described  above  and  the  expression  in  the  denominator  was  solved  using  equation  (45) 
with  MLE  estimates  used  for  the  component  Weibull  parameters. 

Figure  33,  Figure  34,  and  Figure  35  show  the  resulting  cost  verses  T  scatter  plots  for 
the  ARP,  BRP,  and  ORP  respectively.  The  respective  true  cost  functions  are  shown  for 
comparison.  The  figures  suggest  the  empirical  cost  models  very  accurately  reflect  the 
theoretical  results  with  very  little  variability.  Figure  36  shows  the  ORP  optimal  value  of 
T  versus  T  scatter  plot.  Again  the  empirical  model  performed  very  well  compared  to  the 
theoretical  result. 


xxx  Data 
Mean 
"  ”  True 


Figure  33.  ARP  Cost  vs.  T 
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IV.  AMRAAM  Data  Analysis 


Overview 

In  this  chapter  the  AMRAAM  failure  data  is  analyzed  with  the  methods  developed  in 
Chapter  3.  As  of  this  writing,  the  AMRAAM  failure  data  consists  of  n  =  815  total 
observations,  606  of  which  are  censored  resulting  in  74.4%  random  right  censoring.  The 
data  set  consists  of  system  failure  data  only.  Lacking  component  level  data,  the  ARP  was 
the  only  preventive  maintenance  policy  that  could  be  considered  for  analysis. 

The  chapter  begins  with  an  analysis  of  the  empirical  MRL  function.  The  empirical 
survivor  functions  are  developed  and  the  non-parametric,  parametric,  semi-parametric 
MRL  functions  are  examined.  The  chapter  concludes  with  an  application  of  the 
reliability  cost  model  developed  in  Chapter  3. 

MRL  Analysis 

Survivor  Function  Estimates.  Figure  37  shows  the  PEXE,  Smoothed  KME  {h  = 
100)  and  KME  survivor  function  estimates.  All  three  non-parametric  estimation 
techniques  match  each  other  very  closely.  Of  note  is  the  relative  “smoothness”  of  the 
KME  survivor  function.  The  very  small  “steps”  in  the  KME  survivor  function  are  due  to 
the  compactness  of  the  data  for  t  <  1600.  Most  startling,  however,  is  the  obvious 
truncated  nature  of  all  three  of  the  non-parametric  survivor  function  estimates.  The 
survivor  function  estimates  abruptly  end  at  r  =  Z„  =  2070  with  5(Z„)  =  0.3.  A  Weibull 
distribution  was  fitted  to  the  data  with  MLE  parameters  (4  =  0.000616  and  =  0.853 
determined  via  equation  (48).  The  Weibull  parametric  survivor  function  estimate  is 
compared  to  the  KME  estimate  in  Figure  38.  The  non-parametric  estimate  is  a  close 
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match  to  the  parametric  result  for  t  <  Zn,  but  the  lack  of  information  in  the  tail  of  the  non- 
parametric  distribution  is  evident. 

KME  PEXE 


Figure  37.  Non-Parametric  Survivor  Function  Estimates 
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—  KME 
■  ■  Parametric 

Figure  38.  ICME  and  Parametric  Survivor  Function  Estimates 

MRL  Function  Estimates.  Figure  39  shows  the  KME,  semi-parametric  KME,  and 
Weibull  parametric  estimates  of  the  MRL  function.  The  smoothed  KME  and  PEXE 
MRL  estimates  are  time  intensive  to  compute  and  were  not  pursued  due  to  the  already 
relatively  smooth  nature  of  the  KME  MRL  function.  The  figure  illustrates  the  dramatic 
departure  of  the  non-parametric  BCME  MRL  function  from  the  parametric  and  semi- 
parametric  results.  The  decreasing  MRL  nature  of  the  non-parametric  KME  result  is 
most  certainly  exaggerated  since  the  non-parametric  KME  survivor  function  contains  no 
information  in  the  tail  of  the  underlying  distribution.  The  semi-parametric  and 
parametric  results  are  in  close  agreement  and  indicate  a  slightly  increasing  MRL 
function.  However,  these  results  may  also  differ  significantly  from  the  true  underlying 
distribution  since  such  a  large  portion  of  the  survivor  function  is  estimated  beyond  the 
data  as  shown  in  Figure  38. 
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Semi-parametric  KME 
■  ■  Parametric 
-  KME 

Figure  39.  KME,  Semi-Parametric  KME,  and  Parametric  MRL  Functions 

MRL  Statistical  Tests.  Chen  Hollander  and  Langberg’ s  V  test  and  Lim  and  Park’s 
S„‘^  test  were  applied  to  the  data  to  test  for  DMRL  and  NBUE  respectively.  The  V  test 
resulted  in  a  test  statistic  value  of  1.967  and  a  corresponding  p- value  =  0.025.  The  S„^ 
resulted  in  a  test  statistic  value  of  5.855  with  a  corresponding  p-value  «  0.001.  Both 
tests  confirm  the  DMRL  /  NBUE  impression  of  the  KME  MRL  function.  However,  these 
tests  are  unreliable  under  the  sample  size  and  censoring  percentage  conditions  of  the 
AMRAAM  data  set  as  demonstrated  in  Figure  18.  Based  on  the  impression  of  the 
parametric  and  semi-parametric  MRL  functions,  I  would  not  reject  Ho  in  favor  of  DMRL 
or  NBUE  under  these  conditions. 

Preventive  Maintenance  Policies 

The  aging  criterion  must  be  satisfied  before  a  preventive  maintenance  policy  is 
pursued.  The  MRL  analysis  presented  above  does  not  support,  nor  necessarily  reject,  this 
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criterion.  Without  conclusive  evidence  of  aging,  the  pursuit  of  a  preventive  maintenance 
policy  for  the  AMRAAM  based  on  the  current  data  is  not  warranted.  However,  for 
illustration  purposes  the  reliability  cost  model  developed  in  Chapter  3  is  applied  here 
using  the  non-parametric  MRL  function. 

Reliability  Determination.  The  system  reliability  S(Tv)  was  computed  for  values  of 
T  from  500  to  2000  hours  in  200  hour  increments  using  the  KME  survivor  and  MRL 
functions.  Figure  40  shows  the  resulting  reliability  versus  T plot. 


Figure  40.  Reliability  vs.  T 


ARP  Cost.  The  normalized  ARP  cost  given  by  equation  (53)  was  computed  for 
values  of  r  from  500  to  2000  hours  in  200  hour  increments.  Figure  41  shows  the 
resulting  ARP  cost  versus  Tplot.  Note  the  dramatic  increase  in  cost  for  T  <  500  which 
represents  a  desired  reliability  >  0.9  as  shown  in  Figure  40.  Figure  42  shows  the  direct 
relationship  of  the  ARP  cost  with  the  desired  reliability  goal. 
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V.  Summary  and  Conclusions 


Research  Objectives  and  Primary  Results 

The  overall  objective  of  this  research  effort  was  to  formulate  a  preventive 
maintenance  strategy  for  using,  retiring,  or  refurbishing  AMRAAM  missiles  subject  to 
extended  captive  carry  flight  time.  A  preventive  maintenance  policy  is  only  applicable  if 
the  item  in  question  is  aging,  or  deteriorating  with  time.  Therefore,  a  supporting  objective 
of  this  research  is  to  characterize  the  aging  process  of  the  missile  system  through  a  non- 
parametric  analysis  of  its  Mean  Residual  Life  (MRL)  function. 

The  MRL  analysis  of  the  current  AMRAAM  failure  data  did  not  support  the  aging 
criterion.  Without  conclusive  evidence  of  aging,  the  pursuit  of  a  preventive  maintenance 
policy  for  the  AMRAAM  based  on  the  current  data  is  not  warranted.  Instead  a  failure 
replacement  maintenance  policy  whereby  individual  components  are  repaired  /  replaced 
upon  failure  is  recommended. 

Summary  of  Other  Signiticant  Results 

Non-parametric  MRL  Functions.  Three  non-parametric  MRL  function  estimation 
techniques  discussed  in  the  literature  were  examined  via  a  numerical  example.  All  three 
estimation  techniques  differed  significantly  from  the  true  MRL  function  in  that  they 
exhibited  a  greatly  exaggerated  decreasing  MRL.  The  departure  from  the  true  result  was 
attributed  to  the  lack  of  information  in  the  tail  of  the  respective  non-parametric  survivor 
functions.  A  semi-parametric  technique  for  estimating  the  MRL  function  was  developed 
whereby  the  tails  of  the  respective  non-parametric  survivor  functions  were  estimated  with 
a  parametric  result. 
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The  semi-parametric  MRL  functions  showed  dramatic  improvement  over  the  non- 
parametric  estimates. 

MRL  Statistical  Tests.  The  performances  of  several  statistical  tests  for  the 
distribution  classes  of  DMRL  and  NBUE  were  assessed  via  Monte  Carlo  simulation. 
Complete  data  and  censored  data  tests  were  considered.  The  simulations  revealed  the 
behavior  of  the  tests  with  changes  in  sample  size,  amount  of  censoring,  and  the  MRL 
characteristic  of  the  underlying  distribution.  The  simulations  also  showed  that  none  of 
the  three  censored  data  tests  considered  were  reliable  under  the  conditions  of  1)  sample 
size  of  815  with  75%  random  right  censoring,  or  2)  sample  size  of  2500  with  80% 
random  right  censoring. 

Maintenance  Policy  Optimization.  Three  preventive  maintenance  policies 
discussed  in  the  literature  were  considered.  The  traditional  approach  of  optimization 
through  minimizing  the  cost  functions  requires  the  cost  for  a  system  failure  (cf,)  be 
explicitly  known.  However,  Cb  is  often  subjective  and  difficult  to  express  in  monetary 
terms.  A  “reliability  cost  model”  was  developed  whereby  system  reliability  for  each 
policy  was  expressed  as  a  function  of  cost.  This  technique  allows  an  informed  decision 
on  the  “optimal”  policy  based  on  a  subjective  assessment  of  Cf,.  Theoretical  results  were 
presented  and  the  performance  of  the  model  applied  to  empirical  data  was  assessed  via  a 
Monte  Carlo  simulation.  Finally,  the  sensitivities  of  the  respective  policy  costs  with 
changes  in  the  cost  parameters  were  examined.  This  analysis  revealed  the  nature  of  the 
relative  policy  costs  as  summarized  below: 

(1)  The  most  cost  efficient  policy  is  irrespective  of  the  desired  reliability  goal. 
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(2)  The  cost  ratio  cjY^Ci  =  0.75  represents  the  boundary  between  the  ARP  and 
BRP  costs.  The  ratio  c.,/Ec,  >  0.75  favors  the  ARP,  while  cjYjCi  <  0.75  favors 
the  BRP. 

(3)  The  ORP  is  always  at  least  as  cost  efficient  as  either  the  ARP  or  the  BRP.  The 
ORP  is  degenerate  with  the  ARP  for  large  values  of  c/Zc,  and  is  degenerate 
with  the  BRP  when  decreases  to  zero.  The  ORP  offers  the  greatest 
advantage  in  cost  efficiency  over  the  other  policies  when  Cj/Sc,  =  0.75. 

Suggestions  for  Future  Research 

MRL  Statistical  Tests.  The  MRL  statistical  tests  considered  in  this  research 
effort  were  determined  to  be  unreliable  when  the  data  is  subject  to  heavy  random  right 
censoring  (>  75%),  even  when  the  sample  size  is  extremely  large  {n  =  2500).  A  statistical 
test  that  it  is  more  robust  under  these  conditions  is  required.  New  tests  motivated  by  the 
notion  of  virtual  age  could  be  pursued.  Recall  that  0  <  t^it)  <  n  for  all  t  for  an  item  with 
a  NBUE  distribution  and  tv(t2)  ^  tvih)  for  all  ^  0  if  an  item  is 
characterized  with  a  DMRL  function. 

AMRAAM  Preventive  Maintenance.  The  MRL  analysis  of  the  AMRAAM  failure 
data  did  not  support  the  aging  criterion.  However,  there  was  great  disparity  between  the 
non-parametric  and  parametric  results  due  to  the  lack  of  data  in  the  tail  of  the  non- 
parametric  survivor  function.  This  subject  should  be  revisited  as  more  data  becomes 
available  allowing  better  approximations  of  the  system  MRL  function  . 

The  AMRAAM  failure  data  available  for  this  research  did  not  contain  component 
level  information.  It  is  well  known  that  complex  systems  often  exhibit  exponential 
lifetime  distributions  even  though  one  or  more  components  may  be  aging.  Component 
MRL  functions  should  be  developed  as  data  becomes  available.  Preventive  maintenance 
policies  at  the  component  level  could  then  be  pursued  if  warranted. 
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Appendix  A:  MRL  Example  Code 


I.  Preliminaries 

1 .  Set  parameters  and  constants 

Sample  size:  n=30 

Weibull  parameters:  p^i.s  a=i 

Plot  time  parameters:  tmax=3.c  tstart  =0.001  tstep=.05 

Smoothing  parameter:  h=o.i 

Non-parametric  or  “semi-parametric"  results:  switch  =0 

switch  =  0,  non-parametric 

switch  =  1 ,  semi-parametric 

2.  Weibull  MRL  and  survivor  functions  (Leemis  1995:88): 

The  incomplete  Gamma  Function  (Leemis  1995:286): 

I(a,x):  =  — !-•  y“"‘-exp(-y)dy 

r(a)  Jo 

Survivor  function: 

S(t, shape  .scale)  :=exp[-(scale-t)®’’®’’®] 

MRL  function: 

m( t , shape , scale )  :  =  exp[(scale-t)^'^=‘P^]  ^ L  _  f _1_ ^ ^ .^^shape 
scale -shape  \  shape  /  [  [shape 


t  := 0.001,0.02..  tmax 


1 

^ ^ - 

2 

- 1 - 1 - 

S(t,P,a) 

m(t,P  ,a) 

y' 

S(_t,l,a)  0.5 

^t,0.8,a) 

0 

m(t,l,a)  1 

m(t,0.8,a) 

n 

_ 1 _ 1 _ 

0  ^ ^ ^ —  0  ^ - - - 

0  1  2  3  0  1  2  3 

t  t 


Weibull  mean  calculation:  Utme  -^ — i^trae=®-9®3 

p-a  \p/ 

3.  Weibull  random  variate  generation: 

X  ;=rweibull(n,P)-a  U:  =  sort(X)  meanX  :  =  mean(X) 
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4.  Weibull  Parameter  Estmation  (Leemis,  1995:216): 
a)  Create  set  of  ordered  observed  failure  times: 


oft(G) 


i4-l 

out^G  if  cols(G)= 
while  i<n 


if  G  =1 
1 


out  ^ 
1 

i<--i'h  1 


G 


,0 


1 

if  cols(G)=2 


out 


Pn(G) 


b)  shape  parameter  calculation: 
estimation  tolerance:  e  :=o.oi 


x<-ofl(G) 
r<-rows(x)—  1 


— -h  5] 

Pn  ^  U, 


i=l 


i=l 


S  prime^  « 

Pn 


P  new^P  n“: 


i=l 


i=l 


i2 


i=  1 


i=l 


» prime 


while  |Pnew-Pn|>£ 

new 


n 


^  i=i 


^  £  ho) 

i=  1 _ 

£  (-=,.0)'" 

i=l 


S  prime^  2 


P  new^P  n 


£  (\o) 

i=l 


Sn 


i=l 


i=l 


£  (“m 

i=l 


®  prime 
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c)  scale  parameter  calculation: 

P„(Z)=  1.222  b:=P„(Z) 

a„(G,p):=  num^n  if  cols(G)sl  a ^(z.p  jj(Z))  =0.953  a  :=a  n(z,p  „(Z)] 

num«-ZG^*^  if  cols(G)»2 


1 


II.  Complete  Data  Example 

1.  Sample  Survivor  Function  (Leemis,  1995:253): 

s„(t)  -li-i 


“  '  True 


2.  Compute  the  semi-parametric  MRL  numerator  constant: 


C:= 


C<— 0  if  switch=0 
rlOZ 

n 

C<-  S(u,p,a)du  if  switch=l 

Z 

•'  n 

c 
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3.  Sample  MRL  Function  (Guess  and  Proschan,  1985:8-9): 


mjj(t,G)  := 


j-1 

while  j^n 
if  t<G 

J 


k^j-1 


mrl^- ^ 


i=k-h  1 


-hC 


n  -  k 


j«-n-h  1 

jH+1 

mrl 

4.  Smooth  Empirical  MRL  Function  (Kalasekera,  1991): 
Integrated  Epanechnikov  kernal  function: 

W  ]j(u)  :=  0.3354102U-  0.0223607u%  .5  if  |u|  <^5 
1  if  u>^ 

0  if  uK-^fs 

smoothed  MRL  function: 


mnk('t.G.h) 


out<-0  if  T>G 


meanX- 


out<- 


1-1  £  w, 

j=l 


x-Gl 


dx 


-hC 


in-  1 


T-Gl 


W, 


\i=l 


if  T<G 


out 
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5.  Compute  MRL  Plot  array 

’^Lcmpit:=  It^tstart 
i<-0 

I  while  t^tmax 

OUtjo-t 

out  j<-nijjj,(t,Z,h) 

outj  2<-nijj(t,Z) 

outj  3^m(t,P,a) 
t  -h  tstep 
1 

out 


i:=0..rows(MRLj.^pl^)-  1 


(^MRL 

(mrl 

(mrl 


6.  Hollander  and  Proschan  (1975)  V*  Test: 


c(k) 


star 


T/^i3  \  3^ 

±t_4.„.k!L3.„2.k_iL 

2  j2' 

^  n  k 

[[\  3  /  2. 

2  2_ 

n 

T  E 

i=l 

k 

-I-- 

6 


meanX 


test  :=(210nr^-Vgt^j. 


p value  :  =  1  -  pnorm(test  ,0, 1) 


test  =0.807 


pvalue  =  0.21 


7.  Aly  (1990)  Tn  Test: 


j=l 


n  M-J 


j  =  i 


j-  1 


1 


meanX 


£ 

i=  1 
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J. 

n^-[t  -meanx] 

T  jj  :  = - - - - L  Tj^  =0.839  pvalue  1-  pnorm^T  j^,0, 1 

2 


pvalue  =  0.201 


asqr  ^  -meanX 

8.  Ahmad  (1992)  Un  Test 


Upval(G)  := 


Un-0 

for  je  2..n 

j 


U, 


1 


n! 


(n-2)! 


test^(3-n)^-Uj^ 

pval^  1  -  pnorm(  |  test  |  ,0,1) 
pval 


Upval(Z)  =0.037 


Censored  data  Example 

1 .  Generation  of  Random  Censored  Data; 
exponential  censoring  distribution: 


Yl:=rexp(n,0.8) 


censored  data  pairs  function:  rndcensor(Gi,G2)  := 


i^O 

while  i<n  -  1 


if  G1.<G2 
1  1 

if  Gl>G2 
1  1 

i^i4-  1 


U 
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function  to  add  a  leading  zero  to  the  data  set:  addO(G)  ;= 


while  i<n 


Zi.o-G-I.o 

—  i  'f-  1 


the  ordered  random  censored  data  set:  zc  :=add0(csort(rndcensor(x,Yi),0)) 


compute  censor  percentage:  percentcensor  -  percentcensor  =  0.467 


last  observation  is  considered  “observed:” 


Zc  ,  :=  1 

n,  1 


2.  KME  Survivor  Function 


1  if  t^Gj  Q 

if 


i^2 

while  i^n 


if  t<G 


1,0 
k^i 
1 

for  je  1,2.. k~  1 
^ — L  if 


n- j-h  1 


i^n-h  1 


s<— 0  if  t>G 


n,0 


3.  Smoothed  KME  Survivor  Function  (Kulasekera,  1990): 

step  size  at  each  observation  of  SKME: 

step(G,i)  :=  I S  kme(^’^,o)  ^  ^ KME(^’^i.Hl,o) 
PkMe(*^’^,o) 

the  smoothed  survivor  function: 
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sum<— 0 
for  je  l..n 

sum^sum-h  W  |^| 

S<—  1  -  sum 
S 


!■  / 


•step(G,j)  if  G  al 


4.  PEXE  Survivor  Function  (Kim  and  Proschan,  1991): 


oft(G) 


j^l 

while  i<n 


if  G.  j=l 

out. 

J.l 

i<—  i-f  1 


out 

compute  sample  hazard  function: 

1 


rn(G,k)  :=- 


oft(G)^  1-  1 

£  ("-‘’(q+i.o-Ci.o) 

i  =  of,(0)^_ 


compute  PEXE  of  S(t): 
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i^2 


s^exp((-rn(G,l)-t))  if  t<oft(G)j  j, 
Zoft<-oft(G) 


if  Z.ff  <t<G  _ 
oft  „,o 


while  i<n 


if  t<Z 


oft. 


,0 


i-  1 


lamdat^  j]  r„(G,k)-JZoft^^  Zoft._j 

k=  1 

i<— n  -h  1 
i<—  i  -H  1 

s<— exp((-lamdat)) 


s<— 0  if  t>G 


n,0 


5.  Compute  Survivor  Curve  data  array 

Srvr:=  It^tstart 
i^O 

while  t^^tmax 
outj  o^t 

outj  i<-SpEx^Zc,t) 
outi_2^SsKME(Zc,t) 

°“'^i,3^^KME^Zc,t) 

outj  S(t,P,a) 

t  -j-  tstep 
i<—  i  1 
out 


i  :=0..rows(Srvr)  -  1 


—  PEXE 

Smoothed  KME 

KME 

Parametric 


6.  Compute  Semi-parametric  MRL  Function  Numerator  constant 

compute  Weibull  parameter  estimates:  a:=an(zc,pj,(Zc))  a  =  1.075 

b;=Pjj(Zc)  b  =  1.315 


(G,i) 
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C<-0  if  switch=0 


lazc 


n,0 


Zc 


n,0 


S(u,b,a)du  i 


C 


switchssl 
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7.  KME  MRL  Function  (Chen  et.  all,  1983): 


i<-  1 

while  i<n-  1 

S4-S 


KME' 
denom<-  S 


G,  G 


i,0. 


int<—  S- 


for  je  i-h  l,i'h  2..n 
S^S 

int^int-h  S-^G 

denom 
i^n  -i-  1 
i<— i-1-  1 


KME(G’^j,o)  if  Gj- 1,1*1 


G. 


j-i.o, 


mrk- 


SKME(G.G„,o)'(fln,o-‘)  +  f^ 


^KME^G.t) 


if  G  ,  .<t<G  . 

n— 1,0  n,0 


mrl 


8.  Smooth  KME  MRL  Function  (Kulasekera,  1990): 


rG 


I^sKME^G)  ' 

”^sKME(G-t,M-,C)  := 


n,0 


^sKME(G.u)du 


M-  ■-^^sKME(^‘^) 


0 


^"“^sKME^G.t) 

if  switch=0 

out<— 0  if  S^O 


out 


if  switch=l 


out^O  if  t^G 


1 

rt  ^ 

r 

^sKME(G,u)du 

+-C 

\  J 

0  I 

otherwise 


n,0 


out 


1 

rt  \ 

r 

^sKMe(G.u)  du  H-  C 

\  J 

0  / 

otherwise 


out 
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9.  PEXE  MRL  estimate  (Joe  and  Proschan,  1981): 


SpExd^Ci,y)  dy  -I-  C 

mrfc — ^ -  if  t<G  „ 

mrl^O  otherwise 
mrl 


10.  Compute  MRL  Plot  Data  Array 


^^^cnsrd 


t<—  tstart 
i^O 


while  t^tmax 
out  o-t 

,  3^  KME^  Zc ,  t ,  C) 
out.  4^m(t,(3,a) 

t<- 1  -h  tstep 
i^i-h  1 


out 


(^MRL 

(mrl 

(mrl 

(mrl 
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1 1 .  Chen  et.  all  (1983)  DMRL  test  with  censored  data  c(u)  = 


n  -  u 


n  -  u  -h  1 


i=  1 


i- 1  /  i- 1 

Zc. 


-i  n  ^ 


j=i 


S^l 

for  ie  2..n 


7  n  '<'> 

\  j=i 


2  Zc.  ,  1 


i-  1 


int^int+S-(G.  p-G._j  J 


1  I — r  4-Zc.  . 

■=(»  '■ 


(^i,0 


j=l 

M-hat  jj 


int 


APn  /  a 

Vc:  = -  B(i,a)  :=exp - -Zc 


Hhatn 


Hhat^ 


i,o 


h^lf:-n.fB(n.2)  B(n,3)  ^  B(n,4)  ^  B(n,5)  B(n,6)  ^  B(n,8)\ 

72  18  16  45  18  72  / 


xsqr  := — +  ^ 
720 

i=  1 


n-1  n  /B(i,2)  B(i.3)  ^  B(i,4)  ^  B(i,5)  B(i,6)  ^  B(i,8)\ 

\  72  18  16  45  18  72  j 


(n-  i-i-  l)-(n-  i) 


-i-  half 


test 


1 

n^>Vc 

j. 

2 

Tsqr 


pvalue  ;  =  1  -  pnorm(test  ,0, 1)  pvalue  =  0.24  test  =  0.708 


12.  Urn  and  Koh  (1996)  NBUE  Test : 


n  i-  1  1 

Lc:  =  — - —  [~~j  c(j)  j’’  In  c(j)  +1 

lihat„  "-'ll  11 

"  i=l  j=l  j=l 


W 


half  -n* 


1  ^  ^  (Zc„_„)' 


4  2-lihat 


n/ 


2-iihat  n 


•B(n,4) 


n-  1 


r  £ 


osqr 

4  (n-iH-l)-(n-i) 

i=  1 


1 


Zc 


1,0 


4  2-lthat 


Zc 


n/ 


2|ihat 


•B(i,4)  -  half 
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1 


pvalue  1  -  pnorm(test  ,0, 1)  test  =  3.981  pvalue  =3.433 10 


13.  Lim  and  Park  (1993)  NBUE  Test : 


5„:=£ 

i=i  i  j=i  j=i 


Tsqr<-  halfl-  half2 


pval<—  1  -  pnorm(test  ,0, 1) 
pval 


6pval(G,5  ,p)  :=  halfl^-+  ^  - - - (B(i,4,G,p) -  ^-^(1.3,0,^) |  Ba,2,G,p) 

^  ^  6  (n-i+l).(n-i)  [\  3/2 

half2^n|(B(n,4,G,p)---B(n,3,G,p)]-H--B(n,2,G,p) 


5pval  ( Zc ,  5  ,  phat  ^)  =  0.074 
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Appendix  B:  Statistaical  Test  Comparison  Code 


Complete  Data  Tests 


Tn  (Aly),  Un  (Ahmad),  and  V*  (Hollandar  and  Proschan)  Tests  for 
D(I)MR1_/NBUE  with  Complete  Data 


1 .  Set  Constants: 


n  =  100 

sample  size 

Pmin=0.9 

weibull  shape  param  min 

Pma)^1.5 

weibull  shape  param  max 

deltap^O.l 

weibull  shape  param  step  size 

a  =  1.0 

weibull  scale  param 

nreps  =30 

num  of  repititions  at  each  p 

2.  Hollander  and  Proschan  (1975)  V*  Test  subroutine: 


/-ii3  \  3 

4-k  .  ,2  o  2,  n 

2  u2' 

-h  ^  ^ 

- 4-n-k  -h  3-n  -k - 

l\  3  /  2_ 

2  2_ 

Vpval(G,mean) 


E 

i=  1 


star^ 


mean 


test^(210n)°-^-Vgt3j 

pval<—  1  -  pnorm(test  ,0, 1) 
pval 
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3.  Aly  (1990)  Tn  Test  subroutine: 


Tpval(G,mean)  :  = 


mean 


£ 

i=  1 


Uln  1- 


i-  1 


n  M-J 


In  1 


j=l 

n 


j-1 


osqr 


n4E 


j=i 


l  +  ln  1- 


j-1 


1 

2 

asqr  ^  -mean 
pval<- 1  -  pnorm^T  ^ ,  0, 1 
pval 


4.  Ahmad  Un  Test  Subroutine  (1992) 


Upval(G)  := 


Un-O 
for  je  2..n 

j 


U, 


E 


i=  1 


n.(n-l) 


test<— (3-n)  -Uj^ 

pval^  1  -  pnorm(test  ,0, 1) 
pval 


5.  Add  an  initial  zero  to  the  ordered  data  set: 


addzero(G) 


i<—  1 


while  i:<n 

i<— i-h  1 

Z 
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6.  Compute  p-value  comparison  matrix 


testcompare 


p^Pmin 

i^O 

Pcount<™0 


while  p<Pmax-h^^5^ 
2 


while  i<pcount  -h  nreps  -  1 
data  ^  rweibulK  n ,  p )  •  a 
rneanX^-  mean(  data ) 

7u<r~  addzero(  sort  (data ) ) 

outio-p 

outi  Vpval(Z,meanX) 
outi  Tpval(Z,meanX) 
out.  ^<-Upval(Z) 
i<-i-t-l 
Pcount<— i 
P^  P  -{-  deltap 


7.  Plot  th©  rOSUltS!  i  :=0..rows( testcompare)  -  1 


avg 


Pmax-  pmin 
imax - 

deltap 

for  ie  0, l..imax 
jmin^  nreps -i 
jmax— jmin-h  nreps  -  1 

out.  Q<—  mean( submatrix( testcompare , jmin, jmax,  0, 0) ) 
out .  y  <—  mean(  submatrix(  testcompare , jmin, jmax,  1,1)) 
out.  2^  mean( submatrix( testcompare  Jmin, jmax,  2, 2) ) 
out .  2<—  niean(  submatrix(  testcompare , jmin, jmax,  3,3)) 


j  :-0..  rows  (avg)  -  1 


mean 


mean 


Censored  Data  Tests 


Tn  (Lim  and  Park),  Vc  (Chen,  Hollandar  and  Langberg),  and  Lc  (Urn  and  Koh) 
tests  for  NBUE  /  DMRL  with  Censored  Data 


1 .  Set  Constants: 

n=30  sample  size  Pniin=o.9  weibull  shape  param  min 

1=0.5  exp  censor  rate  Pmax=i.5  weibull  shape  param  max 

a^i.o  weibull  scale  param  deita|3=o.2  weibull  shape  param  step  size 

x,=o.5  exp  censor  rate 
nreps  =5  num  of  repititions  at  each  p 

2.  Chen  et.  all  (1983)  DMRL  test  with  censored  data  : 

c(u)  '■=— — ^  B(i,a,G,mean)  :  =  exp — - — G.  | 

n  -  u-h  1  \  mean  '’  / 


97 


int^  G, 


for  ie  2..n 


S-Skme  G.G.o  if  G_,  =1 


int^int+S-(G^-G_,^, 


VpvaI(G,Vc,H)  :=  halfl^n- 


B(n,2,G,n)  B(n,3,G,H)  B(n,4,G,H)  B(n,5,G,H)  B(n,6,G,n)  B(n,8,G,n) 


B(i,2,G.|i)  B(i,3,G,(l)  B(i,4.G,H)  B(i,5,G,H)  B(i,6,G,H)  B(i,8.G,p) 


{n-i-|-l)(n-i) 


'Bqr< - l-halfl-|-half2 

720 


pval<-  1  -  pnorm(test,0, 1 ) 


3.  Lim  and  Park  (1993)  NBUE  Test  with  randomly  censored  data  : 


5pval^G,5j^,n)  := 


f  - B(i,4,G,n)- 

1=  1 


half2^n-  |B(n,4,G,^)  -  ^•B(n,3,G,^)j -h ^-6(11,2,0,^) 
Tsqr<—  halfl  -  half2 


4-B(i,3,G,^)  B(i,2,G,^) 


pval<—  1  -  pnorm(test  ,0, 1) 
pval 
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4.  Urn  and  Koh  (1996)  NBUE  Test  with  randomly  censored  data: 


Lpval(G,|ii)  :  = 


Lc^-* 


half<-n' 


asqr< — -i- 
4 


/  : 


i=  1 

L\j= 

1 

fi- 

G„,o^ 

-1-  . 

2-\il 

2-\i^ 

PI  C(j)  •  In  P  c(j)  +  1 

/  \  1 j=  1 


o.  '' 


Ko-q-i. 


n-  1 

£ 

i=  1 


■B(n,4,G,n) 

1 


(n-  i+  l)-(n-  i) 


\4  2-H; 


2-H 


B(i,4,G,H) 


-  half 


test^ 


2  , 

n  -Lc 

1 

2 


asqr 

pval<—  1  -  pnorm(test  ,0, 1) 
pval 

5.  Create  Censored  Data  Pairs  : 

rndcensor  (G1,G2)  : 


i^O 

while  i^n  -  1 
if  G1<G2 


if  GL>G2 
1  1 


Ui,o^-G2 
U.  J^O 


addzero(G)  ;  = 


i<-  1 

while  i^n 


i<-ii-  1 


i^i-h  1 


U 


data  (shape ) 


rweibulK  n ,  shape )  •  a 
Y<— rexp(n,^) 
temp<—  rndcensor  ( X,  Y) 
Zctemp^  csort  ( temp ,  0) 
Zc  <—  addzero  ( Zctemp ) 
Zc 

n,  1 

Zc 
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6.  Compute  Censor  Proportion:  cnsrproportion  ( G) 


num<— 

n  -  num 
out< - 

n 

out 


7.  Compute  Products  for  Test  Stat  Computation 


prod  (G,  Stop ) 


outo^l 
OUtj<—  1 
OUt2^1 


if  stop^l 

for  jG  1..  stop 

if  G  =1 

J.i 


Cj-C(j) 

outo^Cj-outj, 

2 

OUt^^Cj  -out 
OUt^^- c  j'^-OUt 


dummy  <—j 


1 

2 


out 
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8.  Compute  P-value  Comparison  Matrix 


testcompare 


Pcount  <-  0 

P<-pmin 

j-o 

while  P<pmax-h 


deltap 


while  j<pcount  -h  nreps  -  1 
Z<—  data(P) 

AF^O 
5Fn-0 
for  is  l..n 

prd^prod(Z,i-  1) 


AF^AF-h 


1 


8Fn-8Fn+(2'P'<l,-P'<‘„)'(Z,,»-Z,-,,„) 


AF 
Vc< - 

5F. 


P 

out. 

J.l 

Vpval(Z,Vc,n) 

out. 

J.2 

5pval^Z,5j,,(xj 

out. 

J.3 

cnsrproportion  ( Z) 

1 

P<~P  -i-  deltap 
pcount 


out 


censorlvl  :  =  mean\testcompare 


<3>^ 


censorlvl  =0.34 
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L 


Vc  P-value 


9.  Plot  the  Results: 


i :  -  0..  rows(testcompare )  -  1 


avg 


Bmax-  6  min 
imax - ' - 

deltaP 

for  iG  0,l..imax 
jmin^  nreps'i 
jma?^  jmin-f-  nreps  -  1 

out.  Q<-  mean( submatrix( testcompare , jmin, jmax,  0,0)) 
out.  j <—  mean(  submatrix( testcompare  , jmin, jmax,  1,1)) 
outi  2^  mean(  submatrix( testcompare , jmin, jmax,  2,2)) 
out 


mean 
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Appendix  C:  Reliability  Cost  Model  Code 


Theoretical  Reliability  Cost  Model 


1 .  Component  Distribution  Input  Parameters 

Weibull  scale  params  a:=(o  0.0003034  0.0002716  o.ooo2848  0.0002736)^ 
Weibull  shape  params  p  :=(o  1.2  1.3  1.4  i.s)^ 

2.  Cost  Coefficients 


component  costs:  c:-(o  .3  .25  .2  .15)^ 

setup  repair  cost:  c  s 

system  replacement  cost:  Cp  -Cg  +  Ic 

r  0 

component  replacement  costs:  c  f :  =  s  ^2 

C  g  +  Cj 


3.  Weibull  Survivor  and  MRL  Functions  (from  MRL  example  worksheet) 

4.  System  Survivor  and  MRL  Functions  tmax=450C  t  := 0.00a  100  .  tmax 
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5.  Function  to  determine  the  Reliability  Goal  given  T 

"Virtual"  replacement  period:  t^CT)  :=mQ(0)-  tno(T) 

System  reliability:  r(T)  :=So(Ty(T)) 


6.  Normalized  Age  Replacement  Policy  Cost  Model 

1 


Ca(T):=- 


So(u)du 


0 


7.  Component  renewal  function: 
Weibull  mean  calculation: 
Weibull  variance  calculation: 
Weibull  density  function: 


H(p,a):=— -r  - 

3-a  3 


var(p,a)  := 


1 


(ar 


P  \P/  \P  \P 


Large  t  renewal  function  approximation: 

W  i(t,P,a)  -  jn^ji(p,a) 

H  J 


\ 


2-11 


s-f(s,P,a)  ds 


Component  renewal  function:  w(t,p,a)  :=max((i- S(t,p,a)  Wj(t,p,a))) 

8.  Normalized  Block  Replacement  Policy  Cost  Model 


It  E  r^-wlT.P,.-,) 


Cb(T) 


i=l  p 
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9.  Normalized  Opportunistic  Replacement  Policy  Cost  Model 


Co(T):= 


T<~0 

"step  ^100 

outo^C^CT) 

OUtj^T 

while  t<T 


denom^t  -h 


rT-t 


exp 


Pi 


i=  1 


1 


new<- 

denom 
if  new<outQ 


out^^—  new 
outj^^T 


1+  £ 


i=l  'P 


T<-T-hT 


Step 


out 


10.  Policy  Costs  vs  Reliability  Tmax:=450C  Tstep:=100 


C0ST:= 


i^O 


T<— Tmax 
while  T>0.001 

out  i^Tv(T) 
out.  2^R(T) 
outj  3<-C2^(T) 
out  4^Cb(T) 
out  5^Co(T)^ 

T<— T-  Tstep 
i<-  i  +  1 


du 
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Policy  Cost 


Cost  Analysis 


Tmax:  =  450C 


Tstep  :=10C 


COST:= 


i<— 0 


T<-Tmax 
while  T>0.001 

out.  J<-C^(T) 
out,  2<— Cg(T) 

out  3^Co(T)^ 

Tf-T-  Tstep 
1 


ARP,  BRP,  ORP  vs  T  Tau  vs  T 


—  ARP 
“  BRP 
"  ORP 


Empirical  Reliability  Goal 

1.  Model  Inputs. 

sample  size:  n=50 

Weibull  param.  estimation: 

shape  param  initial  guess:  p  o=i-3 

shape  param  tolerance:  e=.02 
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exponential  censor  rate:  x,=.ooiJ 

smoothed  KME  bandwidth:  h^ioo 

number  of  replications:  nreps  =4 

replacement  period:  Tmax:=100C  Xstep  :=25C  Tmin  :  =  500 

2.  DATA  Set  Generation. 


DATA(a,P) 


X  1 «—  rweibullfn .  P  i )  •  I  ~  | 

N 

X  2<-  rweibull(n ,  P,)  •  | — ] 

N 

X  3<-  rweibull^n ,  p  •  f— | 


X  4«-  rweibullf  n ,  P.)  •  — 

rJ 

X  «—  augment  ^X  ^ ,  X  2) 

X  augment  ^X  ,  X  3^ 

X^<-  augment  (x  ^ ,  X  4) 

Y<— rexp(n,A,) 
for  ie  1..  n 

X.  „^min(/x^.  X^  X^  X^  Y, 

\\  1-1,0  1-1,1  1-1,2  1-1,3  ‘  V 

X.  if  X.  0<Y._, 

Xj  j<-0  otherwise 


for  je  1..4 
if  X„ 


>X. 


i-l,j-l 

^i.2, 3-1-0 


Otherwise 


X<-csort(X,0) 
for  je  1..4 

\2.j3-l-l 

\l-l 

X 


i-l.j-1 
^-1 
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3.  Empirical  System  Survivor  Functions,  (from  MRL  example  worksheet) 

4.  Weibull  Parameter  Estimation  Functions,  (from  MRL  example  worksheet) 

5.  "Semi-Parametric”  System  MRL  Functions 

rinf 

a)  Numerator  Constant  c(a,p,a)  =  s(t,p,a)dt 

b)  KME  MRL  Function:  (from  MRL  example  worksheet) 

c)  Smoothed  KME  MRL  Function:  (from  MRL  example  worksheet) 

d)  PEXE  MRL  Function:  (from  MRL  example  worksheet) 

6.  Function  to  determine  the  Empirical  Reliability  Goal  given  T 

determine  the  MRL  function  estimation  technique:  switch  =  i 
swithch  =  1:  KME 
swithch  =  2:  PEXE 
swithch  =  3:  sKME 

“  out<-S  mjQ^G,T,C)^  if  switch®  1 

out<— Sp£j^G,|X- mpgj^G,T,C)^  if  switch® 2 

switch®3 

out 
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7.  Reliability  vs  T 


Rel 


n  ’ 


for  ke  0..  nreps  -  1 
Z^DATA(a,p) 

p<— mKME(Z,0,C)  if  switchssl 
)i< — m ppyp^ C)  if  switch"2 
switch=3 

i^O 

Tmax 

while  T>Tmin 
periodj^  .^T 

relj^  .^R j^(Z,T,p,C) 

T^T-Tstep 

i^i-1-  1 

npts^i-  1 

for  je  0..npts 

for  IE  0..  nreps  -  1 

out.  .  .  .^period.  . 

j  •(  nreps) -h  1,0  ^  i,j 

out. ,  .<-rel  . 

j-(  nreps) -1-1,1  i,j 


avg 


Tmax-  Tmin 
imax - 

Tstep 

for  ie  0,  L.imax 
jminc—  nreps  i 
jmax— jminn-  nreps  -  1 
out.  ^<~  mean^submatrix|^Rel  , jmin, jmax,  0, 0 

out.  j^mean^submatrix^ReljjJminJmax,  1,  l' 


half^ 


qt(  0.975,  nreps  -  1)  nreps 


nreps 


nreps  -  1 


out.  j- half-stdev  ^submatrix^Relj^,jmin,jmax,l,lj^ 

out.  out.  j  -h  half-stdev  ^submatrix^Relj^,jmin,jmax,  1, 


out 


no 


Theoretical  reliability  at  each  value  of  T  in  the  empirical  model 

Rel:=  for  ie  0.. rows(avg) -  1 
out.^So(T^(avg.^) 
out 


Data 

Mean 

Lower  95%  Cl 
■O"  Upper  95%  Cl 
True 


Empirical  Policy  Costs 
1 .  Model  Inputs. 

sample  size:  n=3O0 

exponential  censor  rate:  x=.ooiJ 

number  of  replications:  nreps  =4 

Weibull  param.  estimation: 

shape  param  initial  guess: 

shape  param  tolerance:  e=m 

replacement  period:  Tmax:  =  200C  Tstep  =500  Tmin:=500 
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2.  Normalized  Non-parametric  ARP  Cost  Model 


CAn(G,T) 


int<-T  if  0<T<G|  ^ 
T-G  if  T>G  „ 

if 


i<~2 
S^l 
a<~  int 
while  a<T 
S^S 


KME 


if  G 


out<^ 


int-int-^S-(G.o-G._j, 
i^i-h  1 

^  KMe(G>  G  p)  if  G  _  j  ^ 
int^int^-S•(^G. 

2_ 

int 


out 


3.  Normalized  Parametric  BRP  Cost  Model 
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4.  Normalized  Parametric  ORP  Cost  Model 


Con(T,ageT,p,a) 


T 

- 

2 

^  step  ^100 

out  <--C^(T) 


OUt^^T 

while  %<T 

I 


denom<—  t  -{- 


T-t 


expi 


i=  1 


1 


new<- 

denom 
if  new<outQ 


out^<—  new 


out^<~  t 


Cp-f- 


i=  1 


T^-T  -h  X 


Step 


out 
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5.  Empirical  Policy  Costs  vs  Reliability 


COSTj, 


for  ke  0..nreps  -  1 
DATA(a,p) 
i<-0 

T<—  Tmax 
for  Je  1..4 

||Z^,<-submatrix(Z,0,n,2*j,2-j-i- 1) 


P  cn.'^P  n(^  c) 
J  ^  ^ 


while  T>Tmin 


block^i^CBn(T,|icn-acn) 

temp-Co„(T,age^.,p^,„,aen) 

opPki^tempo 
taUj^  tempi 

T^T-  Tstep 


npts^i-  1 
for  je  0..npts 
for  iE  0„  nreps  -  1 


out  /  period.  . 

j-(  nreps) -1-1,0  ^  i,j 

out. ,  .  .  ,<-age.  . 

j-(nreps)-{-i,  1  ®  i,j 

out.  .  .  block.  . 

j-(nreps)-hi,2  i,j 

out.,  7^0pp.  . 

j-(  nreps) -I- 1,3  ^^i,j 

out. ,  ,  ,  .  ,<— tau.  . 

j'Cnreps) -1-1,4  i,j 


out 
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avg  :  = 


imax- 


Tmax-  Tmin 


Tstep 
for  ie  0, l..imax 
jmin^  nreps'i 
jma?^  jmin-h  nreps  -  1 

out.  mean ^submatrix^COST ^ , jmin, jmax,  0, 0^ ^ 

out.  mean ^submatrix^COSTj^, jmin, jmax,  1, 

out.  mean ^submatrix^COST , jmin, jmax,  2,2^^ 

out.  mean^submatrix^COST j^,jmin,jmax, 3, 3^^ 

out.  mean^submatrix^COST ^ , jmin, jmax, 4,4^^ 

out 


Theoretical  ARP  cost  at  each  value  of  T  in  the  empirical  model 


Age  := 


for  ie  O..rows(avg)  -  1 

outi^CA(avg._p) 


1  out 


Theoretical  BRP  cost  at  each  value  of  T  in  the  empirical  model 


Block :  “ 


for  ie  O..rows(avg)  -  1 

outj^CB(avgj  o) 


out 


Theoretical  ORP  cost  at  each  value  of  T  in  the  empirical  model 


Opp 


for  ie  O..rows(avg)  -  1 
out^Co(avgi„)^ 


out 


Theoretical  ORP  value  of  x  at  each  value  of  T  in  the  empirical  model 


Tau  := 


for  ie  O..rows(avg)  -  1 


j  :=0..rows(avg)-  1 


Appendix  D:  AMRAAMData 


z,  5, 

z,  Si 

z,  Si 

z,  Si 

Z  Si 

Z  Si 

Z  Si 

Z  Si 

Z  Si 

1 

1 

14.1 

1 

42.8 

1 

84.5 

0 

178.3 

0 

334.7 

0 

555.1 

1 

840.2 

0 

1611.6 

0 

1.2 

0 

14.2 

0 

43.1 

0 

85.5 

0 

178.7 

1 

334.7 

0 

555.8 

0 

844.4 

0 

1621.4 

1 

1.4 

0 

14.6 

1 

43.2 

0 

86.1 

0 

179.9 

0 

336.1 

0 

563.9 

0 

844.5 

0 

1750.5 

0 

1.4 

0 

15.2 

0 

43.2 

0 

86.4 

0 

180 

0 

346.7 

1 

564 

0 

845.8 

0 

1836 

0 

1.5 

0 

15.2 

0 

43.3 

0 

87.5 

0 

180.7 

1 

346.7 

1 

568.3 

0 

846.7 

1 

1854.3 

0 

1.6 

1 

15.4 

0 

43.3 

0 

87.9 

1 

181.5 

1 

348.8 

1 

571.2 

0 

849.8 

0 

1858.2 

0 

1.8 

0 

16.5 

0 

43.5 

0 

88 

1 

183.1 

0 

353.2 

0 

571.6 

0 

854 

1 

1861.8 

0 

1.9 

0 

16.9 

0 

43.5 

0 

88.3 

0 

190.9 

1 

355.7 

0 

574 

1 

861.9 

0 

1882.6 

0 

2.1 

0 

17 

0 

44 

0 

88.7 

0 

191.1 

1 

355.9 

0 

575.5 

0 

869 

0 

1907.8 

0 

2.1 

0 

17.1 

0 

44.4 

1 

90 

0 

193.8 

0 

357.8 

1 

575.6 

1 

872 

0 

1950 

0 

2.4 

1 

17.6 

0 

44.4 

0 

90.1 

0 

194.2 

0 

358.9 

0 

578.2 

1 

876.1 

1 

2011.3 

0 

2.4 

1 

17.9 

0 

44.5 

0 

92 

1 

194.3 

0 

360.5 

1 

578.2 

0 

877.5 

0 

2012.4 

0 

2.4 

1 

18 

0 

44.6 

0 

92.6 

0 

194.4 

0 

361.4 

1 

578.3 

1 

878.8 

0 

2024.8 

0 

2.4 

0 

18.6 

1 

45.1 

0 

93.4 

0 

194.4 

0 

364.4 

0 

578.9 

1 

883.9 

0 

2054.2 

1 

2.4 

0 

18.7 

0 

45.1 

0 

94.1 

1 

195.4 

0 

366 

1 

582.1 

0 

883.9 

0 

2070.1 

0 

2.7 

0 

19.4 

0 

45.6 

1 

95 

0 

195.6 

0 

367.1 

1 

582.2 

1 

884.3 

0 

2.8 

1 

19.5 

1 

45.8 

0 

95.4 

1 

196.2 

1 

368.8 

1 

582.7 

1 

886.3 

0 

2.8 

0 

19.5 

0 

46.3 

0 

96.2 

1 

199.3 

0 

372.6 

0 

585 

0 

888.9 

0 

2.9 

0 

19.5 

0 

46.4 

0 

96.4 

0 

199.6 

0 

374.7 

1 

585.9 

0 

893 

0 

3 

1 

19.6 

0 

46.9 

0 

96.5 

1 

201.5 

0 

375.9 

0 

591.9 

0 

893.8 

0 

3 

0 

19.7 

0 

47.3 

0 

97.8 

1 

203.4 

0 

377.1 

0 

591.9 

0 

894.4 

0 

3.1 

0 

19.8 

1 

47.3 

0 

99.2 

1 

204.5 

1 

377.3 

1 

592.7 

0 

899.9 

1 

3.4 

1 

20 

0 

47.4 

1 

99.7 

0 

204.6 

1 

377.6 

0 

598.2 

0 

905.6 

1 

3.5 

0 

20.1 

0 

47.8 

0 

100.1 

0 

207 

1 

377.9 

0 

598.7 
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